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Preface 



The design of most modern engineering systems entails the consideration of a 
good trade-off between the several targets requirements to be satisfied along 
the system life such as high reliability, low redundancy and low operational 
costs. These aspects are often in conflict with one another, hence a compro- 
mise solution has to be sought. Innovative computing techniques, such as 
genetic algorithms, swarm intelligence, differential evolution, multi-objective 
evolutionary optimization, just to name few, are of great help in founding 
effective and reliable solution for many engineering problems. 

Each chapter of this book attempts to using an innovative computing tech- 
nique in solving a different engineering problem. The contributions of each 
and every one are summarized in the following paragraph. 

In Chapter 1, the authors approach the traveling salesman using the dis- 
crete differential evolution approach. In Chapter 2, the authors solve reliabil- 
ity optimization problems of series-parallel, parallel-series and complicated 
system, using genetic algorithms. In Chapter 3, the authors present exploit 
particle swarm optimization to build fuzzy systems automatically. In Chap- 
ter 4, the authors propose a maintenance optimization of wind turbine sys- 
tems using intelligent prediction tools. In Chapter 5, the authors apply the 
clonal selection algorithm to economic dispatch optimization of electrical en- 
ergy. In Chapter 6, the authors present novel methods for multi-objective 
evolutionary optimization based on a dynamic aggregation of objectives. In 
Chapter 7, the authors solve the problems related to application mapping into 
a network-on-chip platform, using multi-objective evolutionary optimization. 
In Chapter 8, the authors introduce the theory of chaotic optimization meth- 
ods together with some applications. 

The editors are very much grateful to the authors of this volume and 
to the reviewers for their tremendous service by critically reviewing the 
chapters. The editors would like also to thank Prof. Janusz Kacprzyk, the 
editor-in-chief of the Studies in Computational Intelligence Book Series and 
Dr. Thomas Ditzingcr, Springer Verlag, Germany for the editorial assistance 
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and excellent cooperative collaboration to produce this important scientific 
work. We hope that the reader will share our excitement to present this 
volume and will find it useful. 
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Abstract. Combinatorial optimization problems are very commonly seen in 
scientific research and practical applications. Traveling Salesman Problem (TSP) 
is one nonpolynomial-hard combinatorial optimization problem. It can be describe 
as follows: a salesman, who has to visit clients in different cities, wants to find the 
shortest path starting from his home city, visiting every city exactly once and 
ending back at the starting point. There are exact algorithms, such as cutting-plane 
or facet-finding, are very complex and demanding of computing power to solve 
TSPs. There here, however, metaheuristics based on evolutionary algorithms that 
are useful to finding solutions for a much wider range of optimization problems 
including the TSP. Differential Evolution (DE) is a relatively simple evolutionary 
algorithm, which is an effective adaptive approach to global optimization over 
continuous search spaces. Furthermore, DE has emerged as one of the fast, robust, 
and efficient global search heuristics of current interest. DE shares similarities 
with other evolutionary algorithms, it differs significantly in the sense that 
distance and direction information from the current population is used to guide the 

N. Nedjah et al. (Eds.): Innovative Computing Methods, SCI 357, pp. l |^12.| 
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2 J.G. Sauer et al. 

search process. Since its invention, DE has been applied with success on many 
numerical optimization problems outperforming other more popular 
metaheuristics such as the genetic algorithms. Recently, some researchers 
extended with success the application of DE to combinatorial optimization 
problems with discrete decision variables. In this paper, the following discrete DE 
approaches for the TSP are proposed and evaluated: i) DE approach without local 
search, ii) DE with local search based on Lin-Kernighan-Heulsgaun (LKH) 
method, and iii) DE with local search based on Variable Neighborhood Search 
(VNS) and together with LKH method. Numerical study is carried out using the 
TSPLIB of test TSP problems. In this context, the computational results are 
compared with the other results in the recent TSP literature. The obtained results 
show that LKH method is the best method to reach optimal results for TSPLIB 
benchmarks, but for largest problems, the DE+VNS improve the quality of 
obtained results. 

Keywords: Optimization, Traveling salesman problem, Evolutionary Algorithm, 
Differential Evolution, Variable Neighbor Search, Local Search, Lin-Kernighan- 
Heulsgaun. 

1 Introduction 

Combinatorial optimization problems occur in various fields of physics, 
engineering and economics. Many of them are difficult to solve since they are 
nonpolynomial-hard (NP-hard), i.e., there is no known algorithm that finds the 
exact solution with an effort proportional to any power of the problem size. 

The Traveling Salesman Problem (TSP) is a well-known example of a NP-Hard 
combinatorial optimization problem that involves finding the shortest Hamiltonian 
cycle in a complete graph of n nodes (cities). In the TSP, it is given n cities 1, 
2,..., n together with all the pair wise distances djj between cities i andy. The goal 
is to find the shortest tour that visits every city exactly once and in the end returns 
to its starting city. 

TSPs raise important issues because many problems and practical applications 
in science, engineering, and bioinformatics fields, such as vehicle routing 
problems, integrated circuit board chip insertion problems, scheduling problems, 
flexible manufacturing systems, printed circuit board, X-ray crystallography, 
punching sequence problems, routing, job scheduling problems, and phylogenetic 
tree construction can be formulated as TSPs. 

TSP has been extensively studied. A large number of approaches have been 
developed for solving TSPs (see [l]-[8]). In this context, there is a great interest in 
efficient procedures based on heuristics and metaheuristics to solve it. The TSP 
has received considerable attention over the last two decades and various 
approaches are proposed to solve the problem, such as branch-and-bound, cutting 
planes, 2-opt, 3-opt, simulated annealing, artificial neural network, and taboo 
search [1],[2]. Some of these methods are exact algorithms, while the others are 
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near-optimal or approximate algorithms [3]-[5]. Recently, evolutionary algorithm 
approaches are successfully implemented to the TSP [6]-[8]. In other hand, 
evolutionary algorithms perform a typically incomplete search in the space of 
solutions by iteratively creating and evaluating new candidate solutions. 

Differential Evolution (DE) is a relatively new evolutionary algorithm, which is 
an effective adaptive approach to global optimization over continuous search 
spaces. DE as developed by Storn and Price [9] has proven to be a promising 
candidate to solve real- valued optimization problems [10]. DE has emerged as one 
of the fast, robust, and efficient global search heuristics of current interest. The 
computational algorithm of DE is very simple and easy to implement, with only a 
few parameters required to be set by a user. 

Classical DE approaches use the floating-point (real-coded) representation 
randomly generated initial population, differential mutation, probability crossover 
and greedy criterion of search. 

The applications of DE on combinatorial optimization problems are still 
considered limited, but the advantages of DE include a simple structure, speed to 
acquire solutions, and robustness that are sustained in the literature. However, the 
major obstacle of successfully applying a DE algorithm to combinatorial problems 
in the literature is due to its continuous nature [11]. Aiming at the discrete 
problems, novel discrete DE approaches have been proposed in recent literature to 
solve combinatorial optimization problems [12]-[15]. In other hand, in spite of the 
prominent merits, sometimes DES shows the premature convergence and slowing 
down of convergence as the region of global optimum is approached. The 
application of local search in DE is an alternative strategy to improve the 
convergence performance. 

In this paper, the following discrete DE approaches for the TSP are proposed 
and evaluated: i) discrete DE approach without local search, ii) DE with local 
search based on Lin-Kernighan-Heulsgaun (LKH) [16] method, and iii) DE with 
local search based on Variable Neighborhood Search (VNS) [17], [18] and 
together with LKH method. Computational results evaluated in the TSP are based 
on Reinelt's TSPLIB [19]. In this context, those results are compared with the 
other results in the recent TSP literature. The obtained results show that LKH 
method is the best method to reach optimal results for TSPLIB benchmarks, but 
for largest problems, the DE+VNS improve the quality of obtained results. 

In general terms, the contribution of this paper was the application of an 
extension of DE algorithm using discrete variables and local search mechanisms 
to a set of benchmark problems described in TSPLIB and studied the 
incorporation of different local search schemes to improve the performance of 
discrete DE. Simulation results have been presented to compare the performance 
of different schemes. 

The remainder of this paper is organized as follows. In Section 2, a describtion 
of the fundamentals of TSP is provided. Section 3 presents the features of discrete 
DE approach with local search. Section 4 then describes the TSPLIB benchmark 
problems and evaluates the quality of simulation results. Lastly, section 5 presents 
our conclusion and future research. 
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2 Fundamentals of TSP 

The TSP is one of the most extensively studied problems in combinatorial 
optimization. Problems of combinatorial optimization distinguish themselves by 
their well-structured problem description as well as by their huge number of 
possible action alternatives. 

The task of TSP basically consists of finding the shortest tour through a number 
of cities, visiting every city exactly once. The TSP can be formulated as follows. 
Given matrix D=(d.) and the set n of permutations of the integers from 1 to n, 

find a permutation 7r=(7r(l), n(2), ..., 7r(w))eil that minimizes 

n-\ 
Z(lt)= Z^(,-),^(;'+l) + d 7C ( n ),7t(\) (1) 

The interpretation of n, D and n is as follows: n is the number of cities; D is the 
matrix of distances between all pairs of these cities; j = n{i) denotes city j to visit 
at step i. Usually, permutations are called tours, and the pairs (tt(1), tt(2)), ..., {n{i), 
x(i+l)), ..., (7r(»), 7r(l)) are called edges. So, solving the TSP means searching for 
the shortest closed tour in which every city is visited exactly once. In other words, 
our goal is to find an ordering n or tour, of the cities that minimizes the length of 
the tour, given by Eq. (1). 

In this work, we will restrict our attention to the two-dimensional Euclidean 
TSP, which is the special case where the cities are points in the plane and dy is the 
Euclidean distance from city i to city j. 

If the distances satisfy dy = d ti for 1 < i, j < N, this case is the symmetric TSP. 
However, it is possible to discard that last condition and allow the distance from 
city i to city j to be different from the distance between city j and city i. We refer 
to that case as the asymmetric TSP. 

3 Differential Evolution Algorithm 

A. Classical DE algorithm 

DE is a population-based stochastic function to minimize (or maximize) relating 
to evolutionary algorithms, whose simple yet powerful and straightforward 
features make it very attractive for numerical optimization. 

DE combines simple arithmetical operators with the classical operators of 
recombination, mutation and selection to evolve from a randomly generated 
starting population to a final solution. DE uses mutation which is based on the 
distribution of solutions in the current population. In this way, search directions 
and possible step sizes depend on the location of the individuals selected to 
calculate the mutation values [20]. It evolutes generation by generation until the 
termination conditions have been met. 

The different variants of DE are classified using the following notation: 
DE/ afjB/S, where a indicates the method for selecting the parent chromosome that 
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will form the base of the mutated vector, /? indicates the number of difference 
vectors used to perturb the base chromosome, and Vindicates the recombination 
mechanism used to create the offspring population. The bin acronym indicates that 
the recombination is controlled by a series of independent binomial experiments. 

The fundamental idea behind DE is a scheme whereby it generates the trial 
parameter vectors. In each step, the DE mutates vectors by adding weighted, 
random vector differentials to them. If the cost of the trial vector is better than that 
of the target, the target vector is replaced by the trial vector in the next generation. 
The variant implemented here was DE/ randlMbin, which involved the following 
steps and procedures: 

Step I: Initialization of the parameter setup: The user must choose the key 
parameters that control DE, i.e., population size, boundary constraints of 
optimization variables, mutation factor (fj, crossover rate (CR), and the stopping 
criterion (t ma ). 

Step 2: Initialize the initial population of individuals: Initialize the generation's 
counter, t = 0, and also initialize a population of individuals (solution vectors) x(t) 
with random values generated according to a uniform probability distribution in 
the w-dimensional problem space. 

Step 3: Evaluate the objective function value: For each individual, evaluate its 
objective function (fitness) value. 

Step 4: Mutation operation (or differential operation): Mutate individuals 
according to the following equation: 

Zi (t + 1) = x h {t) + f m ■ [x i2 (t) - x i3 (t )] (2) 

where i =l,2,...,Nis the individual's index of population; / is the time (generation); 
Xj(t) = [Xi. (t),X{. (t),...,X{ (t)\ stands for the position of the ;-th individual of 

population of N real-valued M-dimensional vectors; z ; (t) = [z,- (t), zi ? (t),..., z; (t)\ 

stands for the position of the ;-th individual of a mutant vector, f m > is a real 
parameter, called mutation factor, which controls the amplification of the 
difference between two individuals so as to avoid search stagnation. The mutation 
operation randomly select the target vector x (| (/) , with i & j'j .Then, two 

individuals x i2 (t) and x^ if) are randomly selected with j'i * i^ & i^ * i , and the 

difference vector x,, - x h is calculated. 

Step 5: Crossover (recombination) operation: Following the mutation operation, 
crossover is applied in the population. For each mutant vector, z,(f+l), an index 
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rnbr(i) e { 1, 2, • • • , n } is randomly chosen using a uniform distribution, and a trial 
vector, ^■(r+l)=[w ! ff+lj, u t (t+l),...,Uj (t+l)f , is generated via 

IZi . (t +1) if randb(j) < CR or y = rnbrii ), 
x,- . (f ) otherwise, (3) 

where 7=1,2,..., « is the parameter index; Xy(?) stands for the i-th individual of j-th 
real-valued vector; zijit) stands for the j-th individual of j-th real-valued vector of 
a mutant vector, uy(f) stands for the j-th individual of j-th real-valued vector after 
crossover operation; randb(j) is the y'-th evaluation of a uniform random number 
generation with [0, 1]; CR is a crossover rate in the range [0, 1]; 

To decide whether or not the vector u t (t + 1) should be a member of the 
population comprising the next generation, it is compared to the corresponding 
vector Xj( t ). Thus, if/ denotes the objective function under minimization, then 

r, + n-J M '" ( ' + 1) if /(«.■(' + 1 » </(*i(0). (A , 

Xj\t + 1) —< (4 J 

[Xi(t) otherwise, 

Step 6: Update the generation 's counter. t = t+ 1 ; 

Step 7: Verification of the stopping criterion: Loop to Step 3 until a stopping 
criterion is met, usually a maximum number of iterations (generations), t max . 

B. Discrete DE algorithm 

DE algorithms are evolutionary algorithms that have already shown appealing 
features as efficient methods for the optimization of continuous space functions. 
The DE algorithms use a floating-point representation for the solutions in the 
population. 

However, the continuous nature of the algorithm prohibits DE to apply to 
combinatorial optimization problems. To compensate this drawback, Tasgetiren et 
al. [21], [22] presented the Smallest Position Value (SPV) rule, barrowed from the 
random key representation of Bean [23], for the particle swarm optimization 
(PSO) algorithm, which is developed by Kennedy and Eberhard [24], to convert a 
continuous position vector to a job permutation. It has been successfully applied to 
the single machine total weighted tardiness problem and the permutation flow- 
shop sequencing problem. 

The SPV rule can still be used in DE since smallest position value can be 
replaced by smallest parameter value to convert the continuous parameter values 
to a permutation. Details about the SPV rule for DE's solutions representation are 
presented in [21], [22]. 

For example, with a problem have 6 cities, the first thing is generate a random 
solution, like 1, 3, 5, 6, 2, 4. After this, it is made a mutation of the continuous 
values and found a third vector. This vector is changed to a possible solution, e.g. 
1, 2, 4, 5, 6, 3. If this solution has less length than the first solution, so, this will be 
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the new best solution. And the process will keep continuing until the stopping 
criterion is reached. 

C. Discrete DE with Local Search 

The Lin-Kernighan (LK) heuristic [25] is generally considered to be one of the 
most effective methods (state-of-the-art of local search heuristics) for generating 
optimal or near-optimal solutions for the symmetric TSP. 

Domain-specific heuristics, such as 2-opt, 3-opt, and LK, are surprisingly very 
effective for the TSP. The LK algorithm, also referred to as variable-opf, however 
incorporates a limited amount of hill-climbing by searching for a sequence of 
exchanges, some of which may individually increase the tour length, but which 
combine to form a shorter tour. A vast amount has been written about the LK 
algorithm, including much on its efficient implementation. 

Recently, a new and highly effective variant of the LK algorithm has been 
developed by Helsgaun [16]. This scheme employs a number of important 
innovations including sequential 5-opt moves and the use of sensitivity analysis to 
direct the search. 

Other alternative approach is the local search based on VNS [17], [18]. 
Contrary to other heuristics and metaheuristics based on local search methods, 
VNS does not follow a trajectory but explores increasingly distant neighborhoods 
of the current incumbent solution, and jumps form this solution to a new one if 
and only if an improvement has been made. 

In this paper, the following discrete DE approaches for the TSP are validated: 
i) DE approach without local search, ii) DE with local search based on LKH 
(DE+LKH), iii) DE with local search based on VNS (DE+VNS) and iv) DE+VNS 
combined with LKH method (DE+VNS+LKH). The figure 1 show the tests 
realized. 

In DE approaches with local search mentioned (DE approaches with sequential 
hybridization), the DE is used to explore the more promising part of the TSP 
solution space in order to generate "good" initial solutions, which are refined with 
LKH and/or VNS. 
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Fig. 1 Representation of the proposed models. 
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4 Computational Results 

Each optimization method was implemented in DEV-C++ development platform 
with MinGW compiler (ANSI C/C++ compiler) under Windows XP operational 
system. All the programs were tested in Intel Core Due, 1.6 GHz processor with 
1024 MB of random access memory. In each case study, 50 independent runs 
were made for each of the optimization methods involving 50 different initial trial 
solutions for each optimization method. 

In discrete DE, the population size A' was 50 and the stopping criterion t lnax was 
100 generations (5000 evaluations of the objective function). In the LKH and 
VNS, the stopping criterion t max was 100 generations or if is found best solution. 
In the VNS, the population size was divided in 6 small groups and it was applied a 
local search 2-opt. This local search uses the same stopping criterion of other 
methods. 

The results found with the techniques proposed in this work for the TSPLIB's 
symmetrical cases are compared with optimal value for each tested benchmark. In 
this paper, it considers the results of solving the selected four instances in the 
TSPLIB with the number of cities varying between 101 and 7397. In Tables 1 to 7 
are shown the mean errors/deviations of the obtained results in relation to the 
optimal values of objective function z{7t) and also a statistical analysis of results 

for 50 runs. 

Looking the results, it is possible to verify that in all the instances the DE was 
successful in increase the results. All the mean results show a relevant increasing 
in the performance and the time to find the solution was very fast too. 

In the Tables I and II, it is possible to observe the excellent performance of DE 
+ VNS + LKH, DE + LKH, and LKH approaches. In the Table III, it is possible to 
verify that the standard deviation was better than the original result and the time 
running the DE+VNS+LKH was faster than the original. Table IV shows that the 
DE+VNS+LKH obtained 100% of best solution in all runs. And again the result 
obtained was faster than the original. 

Although, the Tables V and VI, the DE+VNS+LKH had an increasing of time, 
but the standard deviation were less than the others. In neither case the mean reach 
100%, but using the DE was possible to verify that the field was get better. 
Finally, the Table VII shows that in big instances the result found using the DE 
was fast, but the standard deviation increase and either using or not the VNS, the 
value was bigger than the original. 

Table 1 Results in terms of minimization of function z(jt) for the EIL101 problem (50 
runs, Optimum Value = 629). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean 
Time (s) 


DE 


2189 


28.73 


2276 


27.63 


2343 


26.84 


36.11 


1.95 


DE + VNS 


1460 


43.08 


1612 


39.01 


1765 


35.63 


68.35 


0.41 


DE+VNS+LKH 


629 


100 


629 


100 


629 


100 


0.00 


0.02 


DE + LKH 


629 


100 


629 


100 


629 


100 


0.00 


<0.01 


LKH 


629 


100 


629 


100 


629 


100 


0.00 


0.02 
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Table 2 Results in terms of minimization of function z(7Z) for the KROA150 problem (50 
runs, Optimum Value = 26524). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean 
Time(s) 


DE 


168126 


15.77 


177449 


14.94 


182049 


14.56 


3008.65 


2.58 


DE +VNS 


107224 


24.73 


119526 


22.19 


127406 


20.81 


3742.95 


0.67 


DE+VNS+LKH 


26524 


100 


26524 


100 


26524 


100 


0.00 


0.06 


DE + LKH 


26524 


100 


26524 


100 


26524 


100 


0.00 


0.04 


LKH 


26524 


100 


26524 


100 


26524 


100 


0.00 


0.04 



Table 3 Results in terms of minimization of function z(7Z) for the LIN318 problem (50 
runs, Optimum Value = 42029). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean 

Time(s) 


DE 


467121 


8.99 


476169 


8.82 


484940 


8.66 


4353.08 


4.97 


DE + VNS 


317948 


13.21 


337404 


12.45 


359152 


11.70 


8242.69 


1.56 


DE+VNS+LKH 


42029 


100 


42066 


0.99 


42143 


0.99 


53.51 


1.48 


DE + LKH 


42029 


100 


42066 


0.99 


42143 


0.99 


52.21 


0.50 


LKH 


42029 


100 


42067 


0.99 


42143 


0.99 


52.91 


0.50 



Table 4 Results in terms of minimization of function z{7t) for the ATT532 problem (50 
runs, Optimum Value = 27686). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean 

Time 

(s) 


DE 


1343428 


2.06 


1372841 


2.01 


1395667 


1.98 


10863.05 


8.44 


DE +VNS 


901324 


3.07 


965280 


2.86 


1013743 


2.73 


24323.70 


2.66 


DE+VNS+LKH 


27686 


100 


27686 


100 


27706 


0.99 


3.95 


3.34 


DE + LKH 


27686 


100 


27686 


100 


27703 


0.99 


2.40 


1.18 


LKH 


27686 


100 


27688 


0.99 


27706 


0.99 


5.71 


1.38 



Table 5 Results in terms of minimization of function z{7t) for the U2152 problem (50 
runs, Optimum Value = 64253). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean 
Time (s) 


DE 


2358326 


2.72 


2387348 


2.69 


2400930 


2.67 


9857.27 


39.15 


DE +VNS 


1865125 


3.44 


1909209 


3.33 


1954569 


3.28 


17958.15 


11.33 


DE+VNS+LKH 


64253 


100 


64281 


0.99 


64337 


0.99 


23.15 


237.38 


DE + LKH 


64253 


100 


64272 


0.99 


64324 


0.99 


22.60 


83.08 


LKH 


64253 


100 


64278 


0.99 


64324 


0.99 


23.57 


86.58 
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Table 6 Results in terms of minimization of function z{7t) for the PCB3038 problem (50 
runs, Optimum Value = 137694). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean 
Time (s) 


DE 


5133870 


2.68 


5170664 


2.66 


5202107 


2.64 


16979.24 


59.10 


DE + VNS 


4251942 


3.23 


4302541 


3.20 


4356353 


3.16 


26434.48 


17.06 


DE+VNS+LKH 


137694 


100 


137704 


0.99 


137753 


0.99 


20.78 


369.34 


DE + LKH 


137694 


100 


137706 


0.99 


137753 


0.99 


23.34 


123.10 


LKH 


137694 


100 


137710 


0.99 


137757 


0.99 


22.55 


142.64 



Table 7 Results in terms of minimization of function z{7t) for the PLA7397 problem (50 
runs, Optimum Value = 23260728). 



Optimization 
method 


Minimum 


% 


Mean 


% 


Maximum 


% 


Standard 
Deviation 


Mean Time 
(s) 


DE 


2693563800 


0.86 


2710405086 


0.85 


2723186572 


0.85 


6252729.91 


193.23 


DE +VNS 


2282677680 


1.01 


2292893726 


1.01 


2305226654 


1.00 


7921297.27 


131.94 


DE+VNS+LKH 


23260728 


100 


23261170 


0.99 


23264711 


0.99 


1327.66 


13636.55 


DE + LKH 


23260728 


100 


23261533 


0.99 


23265152 


0.99 


1346.24 


6771.52 


LKH 


23260728 


100 


23261000 


0.99 


23264711 


0.99 


1294.50 


7441.50 



5 Conclusion 

In this paper, hybrid discrete DE approaches with local search based on VNS 
and/or LKH were described and evaluated the quality of solutions on instances of 
TSPLIB, and experiments showed its validity (see Tables 1-7). 

As it is possible to verify in the Tables I to IV, the result obtained only using 
the DE was not enough to find a final solution for the TSP problem. But, after 
applied the result found for an initial solution in the LKH problem, it is clear that 
obtained an increasing in the performance. 

It is possible to analyze that for small instances the code is not worth to use, 
but in big instances the results obtained were good and the time to run fall down to 
15% of the original time running only LKH. 

In all the tables, all the configuration parameters were the same, but this could 
be modified to increase more the results. One parameter that could be modified is 
the DE method used. In this paper, it was used the method 1 (DE/rand/1/bin), but 
there are more 9 possible methods based on [10], [1 1]. 

Moreover, the computational results of presented hybrid discrete DE 
approaches with VNS and/or LKH are very close to the best-known TSPLIB 
solution values. Effective implementation of these and related neighborhoods in 
discrete DE approaches are topics for further investigation in multi-objective 
TSPs. 
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The objective of this chapter is to develop and solve the reliability opti- 
mization problems of series-parallel, parallel-series and complicated system 
considering the reliability of each component as interval valued number. For 
optimization of system reliability and system cost separately under resource 
constraints, the corresponding problems have been formulated as constrained 
integer/mixed integer programming problems with interval objectives with 
the help of interval arithmetic and interval order relations. Then the prob- 
lems have been converted into unconstrained optimization problems by two 
different penalty function techniques. To solve these problems, two different 
real coded genetic algorithms (GAs) for interval valued fitness function with 
tournament selection, whole arithmetical crossover and non-uniform muta- 
tion for floating point variables, uniform crossover and uniform mutation for 
integer variables and elitism with size one have been developed. To illustrate 
the models, some numerical examples have been solved and the results have 
been compared. As a special case, taking lower and upper bounds of the inter- 
val valued reliabilities of component as same the corresponding problems have 
been solved and the results have been compared with the results available 
in the existing literature. Finally, to study the stability of the proposed GAs 
with respect to the different GA parameters (like, population size, crossover 
and mutation rates), sensitivity analyses have been shown graphically. 



1 Introduction 

While advanced technologies have raised the world to an unprecedented level 
of productivity, our modern society has become more delicate and vulnera- 
ble due to the increasing dependence on modern technological systems that 
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often require complicated operations and highly sophisticated management. 
From any respect, the system reliability is a crucial measure to be consid- 
ered in systems operation and risk management. When designing a highly 
reliable system, there arises an important question as to how to obtain a 
balance between reliability and other resources e.g., cost, volume and weight. 
In the last few decades, several researchers considered reliability optimiza- 
tion problems, like redundancy allocation and cost minimization problems 
as integer nonlinear programming problems (INLPP) and/or mixed-integer 
nonlinear programming problems (MINLPP) with single or several resource 
constraints [1-14]. To solve those problems, different techniques have been 
proposed by the several researchers. In their works, the reliability of each 
component is known and fixed positive number which lies between zero and 
one. However, in real life situations, the reliability of an individual component 
may not be fixed. It may vary due to several reasons. There is no technol- 
ogy by which different components can be produced with exactly identical 
reliabilities. So, the reliability of each component is sensible and it may be 
treated as a positive imprecise number instead of a fixed real number. Stud- 
ies of the system reliability where the component reliabilities are imprecise 
and/or interval valued have already been initiated by some authors [15-19]. To 
tackle the problem with such imprecise numbers, generally stochastic, fuzzy 
and fuzzy- stochastic approaches are applied and the corresponding problems 
are converted to deterministic problems for solving them. In the stochastic 
approach, the parameters are assumed to be random variables with known 
probability distributions. In the fuzzy approach, the parameters, constraints 
and goals are considered as fuzzy sets with known membership functions or 
fuzzy numbers. On the other hand, in the fuzzy-stochastic approach, some 
parameters are viewed as fuzzy sets/fuzzy numbers and others as random 
variables. However, it is a formidable task for a decision maker to specify the 
appropriate membership function for fuzzy approach and probability distri- 
bution for stochastic approach and both for fuzzy -stochastic approach. So, 
to avoid these difficulties for handling the imprecise numbers by different ap- 
proaches, one may use intervals number to represent an imprecise number, as 
this representation is the most significant representation among others. Due 
to this representation, the system reliability would be interval valued. Here, 
we have considered GA-based approaches for solving reliability optimization 
problems with the interval objective. As the objective function of the relia- 
bility optimization is interval valued, to solve this type of problem by the GA 
method, order relations of interval numbers are essential for selection opera- 
tion as well as for finding the best chromosome in each generation. Here we 
consider the definition of order relations developed by Mahato and Bhunia 
[20] in the context of the optimistic and pessimistic decision maker's point of 
view for maximization and minimization problems. 

In this chapter, we have considered the problem of constrained redun- 
dancy allocation in the series system, the hierarchical series-parallel system, 
the complicated or non-parallel-series system and the network reliability 
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system with interval valued reliability components (redundancy allocation 
and network cost minimization). The problems are formulated as non-linear 
constrained integer programming problems and/or mixed integer program- 
ming problems with interval coefficients [21-22] for maximizing the overall 
system reliability and system cost under some resource/budget constraints. 
During the last few years, several techniques were proposed for solving the 
constrained optimization problem with fixed coefficients with the help of GAs 
[23-29]. Among these methods, penalty function techniques are very popular 
in solving the same by GAs [30-32] . This method transforms the constrained 
optimization problem to an unconstrained optimization problem by penal- 
izing the objective function corresponding to the infeasible solution. Hence, 
to solve the constrained optimization problem the problem is converted to 
unconstrained one by two different types of penalty techniques and the re- 
sulting objective function would be interval valued. So, to solve this problem 
we have developed two different GAs for integer variables with the same GA 
operators like tournament selection, uniform crossover for integer variables 
and whole arithmetical crossover for floating point variables, uniform muta- 
tion for integer variables and boundary mutation for floating point variables 
and elitism of size one but different fitness function depending on different 
penalty approaches. These methods have been illustrated with some numeri- 
cal examples and to test the performance of these methods, results have also 
been compared. As a special case considering the lower and upper bounds 
of interval valued reliabilities of components as same, the resulting problem 
becomes identical with the existing problem available in the literature. 



2 Finite Interval Arithmetic 

An interval number is a closed interval denoted by A — [aL,an] and is de- 
fined by A = [aL,au] = {x : a^ < x < clr,x € 5ft} where cll and an are 
the left and right limits respectively and 5ft is the set of all real numbers. 
A can also be expressed in terms of centre and radius as A = (a c ,a w ) = 
{x : a c — a w < x < a c + a w , x G 5ft}, where a c anda„, are the centre and radius 
of the interval A respectively, i.e., a c = (a^ + a^)/2, and a w = (an — ai)j2. 
Actually, every real number can be treated as an interval, such as for all 
x G 5ft, x can be written as an interval [x, x] having zero width. Now we shall 
present the definitions of arithmetical operations like addition, subtraction, 
multiplication, division and integral power of interval numbers [33] and also 
the n-th root as well as the rational powers of interval numbers [34] . 

Definition 1: Let A = [a_L,a_fj] and B = [&l,6,r] be two intervals. Then the 
definitions of addition, scalar multiplication, subtraction, multiplication and 
division of interval numbers are as follows: 
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• Addition: A + B = [a L , a R ] + [b L ,b R ] = [a L + b L ,a R + b R ] 

• Scalar multiplication: For any real number a,aA = a[a,L,a R ] = 
[cxaL,Oia R ] if a > 
aa R , aa£\ if a < 

• Subtraction: A - B = [ol.Or] - [b L ,b R ] = [a L ,a R ] + [-b R ,-b L ] = 
[a-L - b R ,a R - b L ] 

• Multiplication: Ax B = [a,L,a R ] x [bL,b R ] 

= [mm(a L b L ,a L b R , a R b L , a R b R ),max(a L b L , a L b R , a R b Ll a R b R )] 



• Division 

J_ — u. «_1 V , _ 

Lfor, ' b 



4 = A x -k = [aL,a R ] x [^-,M,providedOi [b L ,h 



Definition 2: Let A = [a,L,a R ] be an interval and n be any non-negative 
integer, then 

< [1, 1] if n = 

.„ _ I [a2' a fl] ^ °i — or ^ n i s oc id 
| I a fl' a il if a -R = and n is even 
I [0, max(a2, a R )\ if a^ < < a/j and n(> 0) is even 

Definition 3: The n — th root of an interval A — [aL, a R ] , n being a positive 
integer, is defined as 

{\Z[a L ,a R ] = [ ^oT , ^or\ if a L > or if n is odd 
[0, ^Ofl] if ol < 0, a/j > and n is even 
if a R < and n is even 
where 4> is the empty interval. 

Again, by applying the definitions of power and different roots of an interval, 

p 

we can find any rational power of an interval. For example Ai obtained by 

p i 

defining Ai as {A p ) q . 



3 Order Relation of Interval Numbers 

Further, for arriving at the optimum solution involving interval algebra, we 
need to define the order relation of interval numbers. 

Let A = [aL,a R ] and B — [&l,6r] be two intervals. These two intervals 
may be one of the following three types: 

1. Type-1: Two intervals are disjoint [see Fig.l]. 

2. Type-2: Two intervals are partially overlapping [see Fig. 2]. 

3. Type-3: One of the intervals contains the other one [see Fig. 3]. 
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Fig. 1 Type-1 interval 



Fig. 2 Type-2 intervals 



Fig. 3 Type-3 intervals 

Here we consider the definitions of order relations developed by Mahato and 
Bfmnia [20] in the context of optimistic and pessimistic decision makers' point 
of view. 



3.1 Optimistic Decision- Making 

In optimistic decision-making, decision maker prefers the lowest value for 
minimization problems and highest value for maximization problems ignoring 
the uncertainty. 

Definition 4: Let us define the order relation > omal between the intervals 
A = [aL,au] and B = [&l,6r] then for maximization problems A >omax B 
&a R > b R ,A > omax B <s> A > omax B KA^B. 

According to this definition, the optimistic decision maker accepts A . The 
order relation > omax is reflexive and transitive but not symmetric. 



Definition 5: The order relation < om in between the intervals A = [aL,an] 
and B = [&z,,6r] then for minimization problems A < om in B <=> cil < &l, 
A < min B <=? A <omin B/\A 7^ B. The order relation < om in is not symmetric. 



18 A.K. Bhunia and L. Sahoo 

3.2 Pessimistic Decision- Making 

In pessimistic decision making, the decision maker prefers the highest/lowest 
value under the principle " Less uncertainty is better than more uncertainty" 
for maximization/minimization problems. 

Definition 6: The order relation > pma x between the intervals A — [aL,aj^] = 
(a c ,a w ) and B — [&l,&.r] = (b c ,b w ) , then for maximization problems 

(i)A > pmax B <=? a c > b c for type — 1 and type — 2 intervals, 
(ii)-A > P max B <=> a c > b c A a w < b w for type — 3 intervals 

However, for Type-3 intervals, pessimistic decision cannot be taken when 
a c > b c A a w > b w . In this case, we consider the optimistic decision. 

Definition 7: The order relation < pm i n between the intervals A — [aL,aj^] = 
(a c ,a w ) and B = [&l,6r] = {b c ,b w ) , then for minimization problems 

(i)A < pm in B o- a c < b c for type — land type— 2 intervals, 
(ii) A < pm in B <=> a c <b c A a w < b w for type — 3 intervals 

However, for Type-3 intervals, pessimistic decision cannot be taken when 
a c < b c A a w > b w . In this case, we consider the optimistic decision. 



4 Assumptions and Notations 

Witout loss of generality, tet us assume the following: 

• The component reliabilities are imprecise and interval valued. 

• The failure of any component is independent of that of the other 
components. 

• All redundancy is active redundancy without repair. 

The following notations have been used in the entire paper. 

the number of redundant con 
reliability of j'-th component 



• Xj : the number of redundant components in j-th subsystem 



• Rj(x): l — (l — rj) Xj ,j = 1, 2, ..., q, the reliability of j-th parallel subsystem 

• x: (xi,X2,—,x n ) 

• r jLi r jR'- lower and upper limits of Tj 

• m: number of resource constraints 

• n: number of stages of the system 

• Rjl (x) : lower bound of Rj (x) 



Rjft(x): upper bound of Rj(x) 



• Rj-. the reliability of j-th subsystem, j = q + l,q + 2,- • •, n 
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(X,R): (X1,X2, —,X q , Rg+i, ..., Rn) 

Rs(x 7 R): system reliability 

Rsl{x, R): lower bound of Rs(x, R) 

Rsr(x,R): upper bound of Rs(x,R) 

Ci(x, R): consumption of i-th resource (i = 1, 2, ..., rn) 

C w (x,R): weighted cost 

cf. availability of i-th resource (i = 1,2,..., to) 

Ij,Uj: lower and upper bounds of Xj 

otj,(3j\ lower and upper bounds of Rj, j = q + 1, q + 2, • • •, n 

R*: minimum prescribed reliability in case of cost minimization problem 

psize: population size 

pjcross: probability of crossover or crossover rate 

pjmute: probability of mutation or mutation rate 



5 Constrained Redundancy Optimization Problem for 
different Systems 

5.1 Series System 

It is well known that a series system (ref. Fig. 4) with n independent compo- 
nents must be operating only if all the components are functioning. In order 
to improve the overall reliability of the system; one can use more reliable com- 
ponents. However, the expenditure and more often the technological limits 
may prohibit an adoption of this strategy. An alternative technique is to add 
some redundant components as shown in Fig. 5. The goal of the problem is 
to determine an optimal redundancy allocation so as to maximize the overall 
system reliability under limited resource constraints. These constraints may 
arise out of the size, cost and quantities of the resources. Mathematically, the 
constrained redundancy optimization problem for such a system for interval 
valued of reliability can be formulated as follows: 

Problem-1: Maximize [R S l,Rsr}=II [{1 - (1 - r,-i)*'},{l-(l - r jR ) x '}] 

subject to gi(x) < Cj, i = 1,2,..., to and lj < Xj < Uj , for j = l,2,...,q, 
where r 3 = [r lL ,r jR ] 

This is a constrained nonlinear integer programming problem with interval 
valued objective. 



Fig. 4 Series System 
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Fig. 5 Parallel series system 

5.2 Hierarchical Series-Parallel System 

A reliability system is called a hierarchical series parallel system (HSP) if 
the system can be viewed as a set of subsystems arranged in a series par- 
allel; each subsystem has a similar configuration; subsystems of each sub- 
system have a similar configuration and so on. For example, let us con- 
sider a HSP system (n = 10, m = 2) shown in the Fig. 6. This system 
has a nonlinear and non separable structure and consists of nested par- 
allel and series system. The system reliability of HSP is given by R$ = 
{1 - (1 - [1 - Q 3 (l - RiRa)]Ri)(l - R 5 Re)}(l - QrQsQ^Rw- Mathemat- 
ically; the constrained redundancy optimization problem for this system for 
interval valued reliability can be formulated as follows: 

Problem-2: Maximize [R S l,Rsr] = {1-{1-(1-[Q 3 l,Q3r\(1-[Ril,Rir} 

[R2L,R2r])) [RiL, -R-1r])(1— [R5L, R5r}[R6L,R6r})}(^~[Q7L,Q7r}[QsL,Q8r\ 
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Fig. 6 Hierarchical series-parallel system 
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[Q9L, Q9r\) [Riol,Rior] subject to g^x) < c { 



i = 1,2, 



. . , m and lj < Xj < 



for j — 1, 2, ..., q. This is an INLP with interval valued objective. 



5.3 Complicated System 

When a reliability system can be reduced to series and parallel configura- 
tions, there exist combinations of components which are connected neither in 
a series nor in parallel. Such systems are called complicated or non parallel 
series systems. This system is also called the bridge system. For example, let 
us consider a bridge system (n = 5, m = 3) shown in Fig. 7. This system con- 
sists of five subsystems and three nonlinear and non-separable constraints. 
The overall system reliability R$ is given by the expression as follows: 

Rs = R 5 (l - QiQa)(l - Q2Q4) + Q 5 [l - (1 - RiR*){l - R3R4)} 

Ri = Ri(xi) and Qi = 1 — Ri for all i = 1, 2, 3, 4, 5 . 

Mathematically, the constrained redundancy optimization problem for such 
complicated system for interval valued reliability can be formulated as follows: 

Problem-3: Maximize [R S l,Rsr] = [R 5 l,R5r}(1 - [Qil,Qir][Q3L,Q3r\) 

(1 - [Q2L, Q2r][Qa Li Qm]) + [Qbl, Q 5 r]{1 - (1 - [Ril, Rir][R2l,R2r})(1 - 
[Rsl, R3r][R4L,R4r])} subject to gi(x) < q, i = 1, 2, ..., m and lj < Xj < 
Uj, ioij = l,2,...,q 



1 2 ~\ 

5 

3~1 1 \ 4 U 



Fig. 7 Complicated system 



5.4 k-out-of-n System 

A k — out — of — n system is a n-component system which functions when at 
least k of its components function. This redundant system is sometimes used 
in the place of pure parallel system. It is also referred to as k—out — of—n : G 
system. An n-component series system is a n — out — of — n : G system 
whereas a parallel system with n-components is a 1 — out — of—n : G system. 
When all of the components arc independent and identical, the reliability of 
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" I n\ 
k — out — of — n system can be written as Rs = X) ■ )f 3 (l — r) n ° , where 

j=k \3 J 
r is the component reliability. 

5.5 Reliability Network System 

Let us consider a network with n subsystems. The goal of the redundancy 
allocation problem is to determine the number of redundant components in 
each of q parallel subsystems and reliability levels of (n — q) general subsys- 
tems so as to maximize the overall system reliability subject to the given 
resource constraints and also to minimize the overall system cost subject to 
the given constraint on system reliability. The corresponding problems are 
mixed-integer nonlinear programming problems as follows: 

Problem-4: Maximize Rs{x, R) = f(Ri(xi),R2(x2), ...,R q (x q ),R q+ i, ...,R n ) 
subject to Ci(x,R) < Ci, i = l,....,m, and 1 < l } < Xj < Uj 7 Xj integer, j = 
1, ...,q, < ay < Rj+i < /3j < 1, j = 1, ..., n — q, 

Problem-5: MinimizeCu,(a;, R) subject to Rs(x, R) > R*, where Rs(x, R) = 
/(i?i(a;i), i? 2 (^2), •••, Rq(x q ), R q+ i,..., R„) 



6 GA Based Constrained Handling Technique 

In the application of GA for solving reliability optimization problem with 
interval objective, there arises an important question for handling the con- 
straints relating to the problem. During the past, several methods have been 
proposed to handle the constraints in evolutionary algorithms [30], [32] for 
solving the same problem with fixed objective. These methods can be clas- 
sified into several types, viz. penalty function techniques, methods that pre- 
serve the feasibility of solutions, methods that clearly distinguish between 
feasible and infeasible solutions and hybrid methods. Among these methods, 
penalty function technique is very well known and widely applicable. In this 
technique, the amount of constraint violations is added /subtracted to the ob- 
jective function in different ways. When the objective function is increased/ 
decreased with a penalty term multiplied by so called penalty coefficient 
there arises a difficulty to select the initial value and upgrading strategy for 
the penalty coefficient. To overcome this difficulty, Deb [30] proposed a GA 
based Parameter Free Penalty (PFP) technique. In this technique, the worst 
fitness value of GA for feasible solutions is considered as the fitness value of 
infeasible solution without multiplying the penalty coefficient i.e., the fitness 
function values of infeasible solutions are independent of the objective func- 
tion value for the same solution. Therefore, according to the PFP technique, 
the converted problem of problem (1-3) is as follows: 
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Maximize [Rsl(x),Rsr(x)] = [Rsl(x),Rsr(x)] 

m m 

- [y^max(0,g^(x) - a), ^max(0,gt(x) - c*)] + 0(x) (1) 

i=l j=l 

where 6{x) = | J^^ ^^ + min [ RsL , RsR ] if x £ S 
and 5 = {x : gi{x) < Cj, i — 1, 2, ..., m and Z < a; < it} 



Here vam[RsL,RsR\ is the value of interval valued objective function of the 
worst feasible solution in the population. Alternatively, the problem may be 
solved with another fitness function by penalizing a large positive number (say 
M which can be written in the interval form as[M, M] ) [18]. This penalty 
function method is known as Big-M penalty and its form is as follows: 

Maximize [R S l(x), R S r(x)} = [Rsl(x), R S r(x)} + 0(x) (2) 

where ^) = { [0 '° ]ifxe ^ 

wnere vyx) | -[R SL ( X ), R SR ( X )] + [_ M , -M] ifx i S 

and S — {x : gi{x) < Ci,i = 1, 2, ..., m and Z < x < u} 



The above problems (1) and (2) are nonlinear unconstrained integer pro- 
gramming problem with interval coefficients. Also, according to the PFP 
technique, the converted problem of problcm-4 is as follows: 

Maximize Rs(x 7 R) — Rs(x, R) 

m 

- Y^ [max(0, Ci(x, R) - c*), max(0, d(x, R) - c t )} + d(x, R) (3) 

i=\ 

0,0] , if(x,R)eX 



where 6(x, R) - ^ _ Rs ^ R) + min ^^ R ^ RjgR ^ R)] lf (Xj R) $ x 

and X = {(x, R) : Ci{x, R) < Ci, i — 1, ..., m and l<x<u,a<R<(3} 

Here min [Rsl(x, i?), Rsr(x, R)] is the value of the interval valued objective 
function of the worst feasible solution in the population. 

Alternatively, the problem may also be solved with another fitness function 
by penalizing a large positive number. The converted form is as follows: 

Maximize Rs(x, R) = R s (x, R) + 9{x, R) (4) 

wherefl(xi?)-I [0 '° ] if{x,R)eX 

wnere u(x, it) - ^ _ Rs ^ R ) + [_ M) _ M ] lf ^ R) ^ x 

and X — {(x,R) : Ci(x,R) < Ci, i = 1, ..., m and I < x <u, a < R < (3} 
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Similarly, for Problem-5, the converted problem is as follows: 
Minimize C w (x, R) — C w (x, R) 

m 

+ J2 [max(0, -R SL (x, R) + R*)} + B(x, R) (5) 



[0,0] if(x,R)€X 

-C w (x,R)+max{C w (x,R)} if(x,R)<£X 

and X = {(x, R) : -R S l{x, R) + R* < 0, % = 1, 2, ..., m and I < x < u, a < R < f3} 



where 6(x,R) 



Here max {C w (x, R)} is the value of the interval valued objective function of 
the worst feasible solution in the population. Alternatively, the problem may 
also be solved with another fitness function by penalizing a large positive 
number. The converted problem is of the form 

Minimize C w (x,R) = C w {x,R) + 6{x, R) (6) 

where 0(x R) = I [ °' ° ] if(x,R)<=X 

wnevev[x,n) < y_ Cw ^ R) + M if( x ,R)<£X 

and 

X = {(x,R) : -Rs L (x,R) + R* < 0, i = l,...,m and I < x < u,a < R < /?} 

The above problems (1-2) are non-linear unconstrained integer programming 
problem with interval coefficients whereas problems (3-6) are non-linear un- 
constrained mixed integer programming problem with interval coefficients. 

7 Genetic Algorithm 

Genetic Algorithm is a well-known stochastic method of global optimization 
based on the evolutionary theory of Darwin: ' The survival of the fittest' and 
natural genetics (Goldberg [23]). It has successfully been applied in different 
real world application problems. This algorithm starts with an initial popu- 
lation of chromosomes. These populations will be improved from generation 
to generation by an artificial evolution process. During each generation, each 
chromosome in the entire population is evaluated using the measure of fit- 
ness and the population of the next generation is created through different 
genetic operators. This algorithm can be implemented easily with the help 
of computer programming. In particular, it is very useful for solving com- 
plicated optimization problems which cannot be solved easily by analytical 
/direct /gradient based mathematical techniques. 
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For implementing the GA in solving the optimization problems, the fol- 
lowing basic components arc to be considered. 

• GA Parameters 

• Chromosome representation 

• Initialization of population 

• Evaluation of fitness function 

• Selection process 

• Genetic operators (crossover, mutation and elitism) 

• Termination criteria 

Initially, the chromosomes/individuals are generated randomly. In this work, 
each chromosome/individual has n components/genes of which first q genes 
are relating to integer variables whereas the last (n—q) are relating to floating 
point variable. These chromosomes/individuals compete with each other with 
their fitness values. Here, the transformed unconstrained objective function 
due to Big-M and PFP penalty are considered as the fitness function. In the 
proposed GA, the well-known tournament selection process is employed as 
the selection operator. The primary objective of this process is to emphasize 
the above average solutions and eliminate the below average solutions from 
the population for the next generation under the well-known evolutionary 
principle " Survival of the fittest" . This selection procedure is based on the 
following assumptions: 

1. When both the chromosomes / individuals are feasible then the one with 
better fitness value is selected. 

2. When one chromosome/individual is feasible and another is infcasible then 
the feasible one is selected. 

3. When both the chromosomes/individuals are infeasible with unequal con- 
straint violations, then the chromosome with less constraint violation is 
selected. 

4. When both the chromosomes/individuals are infeasible with equal con- 
straint violations, then any one chromosome/individual is selected. 

After the selection process, new offspring will be created through crossover 
and mutation processes. In this work, we have used uniform crossover and 
uniform mutation in the genes corresponding to the integer variables, whole 
arithmetical crossover and boundary mutation for the last gene of the 
chromosome. 

The computational steps of crossover are as follows: 

Step-1: Find the integral value of the product of p-cross and psize 

and store it in N . 
Step-2: Select two chromosomes Vk and m randomly from the population. 
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Step-3: For first q genes, compute the components Xkj and Xij 

(j = 1, 2, ..., q) of two offspring by either Xkj — Xkj — g and 

Xij x^j -\- g li Xkj s> Xij or, Xkj — Xj^j -r g ano Xij — Xij g, 

where g is a random integer number between and \xkj — %ij\, 

j = 1 , 2, ..., g and for the last gene, compute the last components 

x' fe and x\a of two offspring will be created by a;jL- = cxkj + 

(f — c)xij and a;^ = (f — cjxkj + cxij where c is a random number 

between and f . 

Step-4: Repeat step-2 and step-3 for ^ times. 

The computational steps of mutation are as follows: 

Step-1: Find the integral value of the product of pjnute and psize 

and store it in N. 
Step-2: Select a chromosome Vi randomly from the population. 
Step-3: Select a particular gene Vik (k = 1, 2, ..., q) of chromosome Vi 

for mutation and domain of Vik is [lit , Uik] . 
Step-4: Create new gene v' ik corresponding to the selected gene Vik 

by mutation process as follows: 

For k = 1,2,..., q 

Vik + A(uik — Vik), if random digit is 
Vik - A(vik -lik), if random digit is 1 

A(y) returns a value in the range [0, y] , is a random integer 

between [0,y]. 

_,,, . , (lik if a random digit is 0. 
Otherwise v' = < . t , , 5 v . , 

I Uik it a random digit is 1. 

Step-5: Repeat Step-2 to Step-4 for N times. 

Sometimes, in any generation, there is a chance that the best chromosome 
may be lost when a new population is created by crossover and mutation 
operations. To remove this situation the worst individual/chromosome is re- 
placed by the best individual/chromosome in the current generation. This 
process is called elitism. The different steps of this algorithm are described 
as follows: 



7.1 Algorithm 

Step 1: Initialize the parameters of genetic algorithm, bounds of variables 
and different parameters of the problem. 

Step 2: t = [t represents the number of current generation]. 

Step 3: Initialize the chromosome of the populationP(t)[ P(t) represents 
the population at t — th generation] . 

Step 4: Evaluate the fitness function of each chromosome of Pit) consider- 
ing any one of the objective function from (1-6) as fitness function. 

Step 5: Find the best chromosome from the population Pit). 

Step 6: t is increased by unity. 
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Step 7: If the termination criterion is satisfied go to step- 14, otherwise, 

go to next step. 
Step 8: Select the population P(t) from the populatiomP(t — 1) of earlier 

generation by tournament selection process. 
Step 9: Alter the population P(t) by crossover, mutation and elitism 

process. 
Step 10: Evaluate the fitness function value of each chromosome of P(t). 
Step 11: Find the best chromosome from P(t) . 
Step 12: Compare the best chromosome of Pit) and P{t — 1) and store 

better one. 
Step 13: Go to step-6. 
Step 14: Print the last found best chromosome (which is the solution of 

the optimization problem). 
Step 15: End. 



8 Numerical Example 

To illustrate the proposed GAs (viz. PFP-GA and Big-M-GA) for solving 
earlier mentioned optimization problems with interval valued reliabilities of 
components, we have solved nine numerical examples. It is to be noted that 
for solving the said problem with fixed valued reliabilities of components, the 
reliability of each component is taken as interval with the same lower and 
upper bounds of interval. For each example, 20 independent runs have been 
performed by both the GAs, of which the following measurements have been 
collected to compare the performances of PFP-GA and Big-M-GA. 

1. Best found system reliability 

2. Average generations 

3. Average CPU times 

The proposed Genetic Algorithms are coded in C programming language 
and run in Linux environment. The computational work has been done on 
the PC which has Intel core-2 duo processor with 2 GHz. In this computation, 
different population size has been taken for different problems. However, the 
crossover and mutation rates are taken as 0.95 and 0.15 respectively. 

Example 1: (related to Problem-1): 

Maximize [R S l,R S r] = II [{1 ~ (1 " r jL ) Xi }, {1 - (1 - r jR ) x '}] subject to: 

£ Pjtf -p<o, 

J=l 

Ec j fe+exp(^)]-C<0, 
i=i 

5 

J2 WjXjexpi^)- W < 0, 
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The values of different parameters along with the interval valued reliabilities 
of Example- 1 are given in Table 1. 



- (1 - [Q 3 l,Q 3 r}(1 

[R5L, R5R][R6L, i?6fl])}(l 



Example 2: (related to Problem-2) 
Maximize [Rsl,Rsr] = {1 - (1 

[RiL,Rir][R2L,R2r]))[R4L,R4r]} (1 — 

[QtL, Q7r}[QsL,Qsr}[Q9L,Q9r\)[RiOL,Rwr] 

subject to 

Ciexp(^)x 2 + c 2 exp(^-) + c 3 x 4 

c 7 a;|exp(^) - 120 < 0, 

wix\x 2 + w 2 exp(^f±) + w 3 x 5 exp(^-) + w 4 x 7 x| + w 5 [x 9 + exp(^>)] 
u! 6 £ 2 exp(2f )-130<0, 



c 4 [ai 5 + exp(^p-)] + c 5 x§x 7 



c 6 x 8 



Table 1 Parameters in Example 1 

3 [rjL,r jR ] pj P Cj C w, W 

1 [0.76, 0.83] 1 7 7 

2 [0.82,0.87] 2 110 7 175 8 200 

3 [0.88, 0.93] 3 5 8 

4 [0.61,0.67] 4 9 6 

5 [0.70, 0.80] 2 4 9 



The values of different parameters along with the interval valued reliabili- 
ties of Example-2 are given in Table 2. 

Table 2 Parameters in Example 2 



3 


[rjL,r jR \ 


Cj 


Wj 


h 


Uj 


1 


[0.80, 0.84] 


8 


16 


l 


4 


2 


[0.87, 0.90] 


4 


6 


l 


5 


3 


[0.89, 0.93] 


2 


7 


l 


6 


4 


[0.84, 0.86] 


2 


12 


l 


7 


5 


[0.88, 0.90] 


1 


7 


l 


5 


6 


[0.90, 0.95] 


6 


1 


l 


5 


7 


[0.80, 0.85] 


2 


9 


l 


3 


8 


[0.90, 0.95] 


8 


— 


l 


3 


9 


[0.80, 0.83] 


— 


— 


l 


4 


10 


[0.88, 0.92] 


- 


- 


l 


6 



Example 3: (related to Problem 2) 

Maximize 

[Rsl.Rsr] = {1 - (1 - 

[RiL,Rir][R2L,R2r]))[R4L,R4r]}(^ — 

[QtL, Q7r][QsL, Qsr][Q9L, Q9r])[RwL, RiQr] 



(1 " [Q3L,Q 3R }(1 - 
[R5L,R5r}[R6L,R6r})}(^ - 
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subject to 

Ciexp(^ L )x 2 + c 2 exp(^.) + c 3 x 4 + c 4 [x 5 + exp(^)] + c 5 x§x 7 + c 6 x 8 + 

c 7 x;jexp(^)-120< 0, 

wixjx 2 + w 2 exp( 2 ^ 1 ) + w 3 x 5 exp(^-) + w 4 x 7 x| + w 5 [x 9 + exp(^)] + 

w 6 x 2 exp(2f )-130<0, 

The values of different parameters along with the interval valued reliabilities 
of Example-3 are given in Table 3. 

Table 3 Parameters in Example -3 



3 


[rjL,r ]R ] 


Cj 


Wj 


h 


Uj 


1 


[0.83,0.83] 


8 


16 


1 


4 


2 


[0.89, 0.89] 


4 


6 


1 


5 


3 


[0.92, 0.92] 


2 


7 


1 


6 


4 


[0.85, 0.85] 


2 


12 


1 


7 


5 


[0.89, 0.89] 


1 


7 


1 


5 


6 


[0.93, 0.93] 


6 


1 


1 


5 


7 


[0.83, 0.83] 


2 


9 


1 


3 


8 


[0.94, 0.94] 


8 


- 


1 


3 


9 


[0.82, 0.82] 


- 


- 


1 


4 


10 


[0.91,0.91] 


- 


- 


1 


6 



Example 4: (related to Problem-3) 

Maximize [Rsl,Rsr] = [R$l, RsrIO- - [Qil, Qir][Q3L, Q3r])(1- [Q2L, Q211] 

[QaLiQar]) + [Q5L,Q5R.]{^ — (1 — [RlL, Rir][R2L, R2r])(1 ~ [R3L,R3R.] 

[R± l ,R A r])} subject to: 

10exp(f-)x 2 + 20x 3 + 3x1 + 8x 5 - 200 < 0, 

10exp(f ) + 4exp(x 2 ) + 2x1 + %4 + ex P(f )] + 7cx p(t) " 310 ^ °> 

12[x| + cxp(x 2 )] + 5x 3 exp(^) + 3xix 



1*4 



4 Jl 1 • '--"•*' v 4 

2xf - 520 < 0, 



(1,1,1,1,1) < (xi,x 2 ,x 3 ,x 4 ,x 5 ) < (6,3,5,6,6), 

where 

fli(si) 

= {[0.78,0.82], [0.83,0.88], [0.89,0.91], [0.915,0.935], [0.94,0. 
i? 2 (x 2 ) = l-(l- [0.73,0.77])^; 

Tit \ ^{Xs + l 

i? 4 (x 4 ) = 1 - (1 - [0.68,0.72])^; 
R^W) = 1 - (1- [0.83, 0.86] ) X6 ; 



], [0.965,0.985]}; 



([0.87, 0.89]) fc ([0.11, 0.13}) X3+1 - k ; 



Example 5: (related to Problem-3) 

Maximize [Rsl,Rsr] = [R5L,R5r\0--[Qil,Qir][Q3l,Q3r])0--[Q2l,Q2r] 

[QaLiQar]) + [Q5L,Q5r]{1 — (1 — [RlL, Rir][R2L, R2r]){1 ~ [R3L,R3R.] 

[R4l,R4r])} subject to 
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10 cxj>(^-)x 2 + 20a;3 + ix\ + 8x 5 - 200 < 0, 

10exp(^) + 4cxp(a; 2 ) + 2x% + 6[xj + exp(^)] + 7cxp(^) - 310 < 0, 

12[x| + exp(x 2 )] + 5x 3 exp(^) + 3x lX j + 2x\ - 520 < 0, 

(1,1,1,1,1) < (x 1 ,x 2 ,x 3 ,X4,x 5 ) < (6,3,5,6,6), 

where 

Ri(xi) 

= {[0.8, 0.8], [0.85, 0.85], [0.9, 0.9], [0.925, 0.925], [0.95, 0.95], [0.975, 0.975]}; 

R 2 (X2) = 1-(1- [0.75,0.75])"*; 

^3(^3) = 3: E 1 ( X3 ^ 1 ) ([0-88, 0.88]) fe ([0.12,0.12])-3+i-fe ; 

i?4(^4) = l~ 2 (l- [0.7,0.7])**; 
R 5 (x 5 ) = l-(1- [0.85,0.85])**; 

The examples 1, 2, 3, 4 and 5 have been solved by two different methods 
PFP-GA and Big-M-GA and the results have been shown in Table 4. 

Table 4 Numerical results for Example 1-5 

Method Exam Popul x Best found system Average Average 

-pie -ation reliability CPU Genera 

size Rs seconds -tion 

(3.2.2.3.3) [0.860808, 0.930985] 0.0001 12.10 
(1, 2, 2, 5, 4, 4, 2, 2, 1, 5) [0.999909, 0.999987] 0.0105 17.55 
(1, 2, 2, 5, 4, 4, 2, 2, 1, 5) [0.999975, 0.999975] 0.0100 17.55 

(5.1.2.4.4) [0.991225, 0.999872] 0.0200 11.20 

(3.2.4.4.2) [0.999382, 0.999382] 0.0100 12.40 

(3.2.2.3.3) [0.860808, 0.930985] 0.0001 12.80 
(1, 2, 2, 5, 4, 4, 2, 2, 1, 5) [0.999909, 0.999987] 0.0110 17.75 
(1, 2, 2, 5, 4, 4, 2, 2, 1, 5) [0.999975, 0.999975] 0.0100 17.75 

(5.1.2.4.4) [0.991225, 0.999872] 0.0200 10.90 
(3,2,4,4,2) [0.999382, 0.999382] 0.0100 12.55 



Example 6: (related to the Problem-4) 

Maximize R s (x, R) = R1R2 + Q2R3R4 + Q1R2R3R4 + R1Q2Q3R4R5 
+Q1R2R3Q4R5 subject to: 



PFP 


1 


50 


-GA 


2 


100 




3 


100 




4 


200 




5 


100 


Big-M 


1 


50 


-GA 


2 


100 




3 


100 




4 


200 




5 


100 



Ci(x) = X1X2 + 2.2a; 2 a;3 + 1.5x2X4 + 2 expf £0i_) < 28, 

C 2 (x) =xi+0.Lr 2 +2:E3 + X4 + 5cxp( 1 M!_) < 25, 

C 3 (x) = x\ + (x 2 -2) 3 + 1.5x 3 +x 4 + 0.6 cxpl^jMf^) < 21, 
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where 1 < Xi < 6, and are integers, i = 1,2,3,4, 0.50 < R5 < 0.99, 
and Ri = Ri(xi) = 1 - (1 - n) Xi ,i = 1,2,3,4, Q i = l- R u i = 1,...,5 
n = [0.69, 0.72] ,r 2 = [0.83, 0.86] ,r 3 = [0.73, 0.76] ,r 4 = [0.79,0.81] 

Example 7: (related to the Problem-4) 

Maximize R s (x, R) = R X R 2 + Q2R3R4 + QiR 2 R 3 Ri + .R1Q2Q3.R4.R5 

+Q1R2R3Q4R5 subject to: 

Ci(x) = x x x 2 + 2.2x 2 x 3 + 1.5x2X4 + 2 exp (^f^) < 28, 

C 2 (x) =xi+0.1x 2 +2a;3 + X4 + 5cxp( 1 %2i_) < 25, 

C 3 (x) = x\ + (x 2 - 2) 3 + 1.5x 3 +x 4 + 0.6 cxp(^) < 21, 

where 1 < Xi < 6, and are integers, i = 1,2,3,4, 0.50 < R5 < 0.99, 
and Ri = R l {x i ) = 1- (l-n) Xi ,i= 1,2,3,4, Q i = l-R u i = 1,...,5 
n = [0.70, 0.70] ,r 2 = [0.85, 0.85] ,r 3 = [0.75, 0.75] ,r 4 = [0.80,0.80] 

The examples 6 and 7 have been solved by two different methods PFP-GA 
and Big-M-GA and the results have been shown in Table 5. 

Table 5 Numerical results for Examples 6-7 



Method Example 


i Popu- 
lation 
size 


(x, R) Best found system Average 
reliability CPU 
Rs seconds 


PFP 6 

-GA 7 


150 
150 


(2,3,1,2,0.9900) [0.958412,0.997223] 0.2705 
(2,1,6,5,0.9396) [0.999927,0.999927] 0.2655 


Big-M 6 
-GA 7 


150 
150 


(2,3,1,2,0.9900) [0.958412,0.997223] 0.2700 
(2,1,6,5,0.9396) [0.999927,0.999927] 0.2590 



Example 8: (related to the Problem-5) 

Minimize C w (x, R) = 0.3Ci(xi) + 0.5C 2 (a; 2 ) + 0.2C 3 (x 3 ) subject to: 
Rs(x,R) > [0.999,0.999], where 1 < Xi < 6, and are integers, i = 
1,2,3,4, 0.50 < R 5 < 0.99, and R s (x,R) , C t {i = 1,2,3) are defined in 
Example 6. 

Example 9: (related to the Problem-5) 

Minimize C w (x,R) =0.3Ci(a;i) + 0.5C 2 (x 2 ) + 0.2C 3 (x 3 ) subject to: 
Rs(x,R) > [0.999,0.999], where 1 < Xi < 6, and are integers, i = 
1,2,3,4, 0.50 < R 5 < 0.99, and R s (x,R) , C,(i = 1,2,3) are defined in 
Example 7. 

The examples 8 and 9 have been solved by two different methods PFP-GA 
and Big-M-GA and the results have been shown in Table 6. 
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Table 6 Numerical results for Example 8-9 



Method Exam Popu- {x, R) 

-pie lation 



Best found Best found system Average 
system cost reliability CPU 

C w Rs seconds 



PFP 8 150 (6,4,2,1,0.8601)33.03866 

-GA 9 150 (2,1,4,4,0.5) 17.97505 



[0.997290, 0.999885] 0.3675 
[0.999081, 0.999081] 0.3525 



Big-M 8 150 (6,4,2,1,0.8601)33.03866 

-GA 9 150 (2,1,4,4,0.5) 17.97505 



[0.997290, 0.999885] 0.3010 
[0.999081, 0.999081] 0.2815 



Table 7 Comparison of results of Ha and Kuo [1] and the proposed methods. 



Example 



System Average 
Reliability CPU 
R s seconds 



Ha and Kuo [1] 4(E2) 

PFP-GA(this work) 3 

Big-M-GA(this work) 3 



(1,1,3,4,2,1,1,3,1,4) 0.999876 

(1, 2, 2, 5, 4, 4, 2, 2, 1, 5) 0.999975 0.0100 

(1, 2, 2, 5, 4, 4, 2, 2, 1, 5) 0.999975 0.0100 



Ha and Kuo [1] 4(E1) (1,3,4,3,3) 

PFP-GA(this work) 5 (3, 2, 4, 4, 2) 

Big-M-GA(this work) 5 (3,2,4,4,2) 



0.999373 

0.999382 0.0100 
0.999382 0.0100 



Table 8 Comparison of results of Sun et al. [9] and the proposed methods. 





Example 


: (x,R) 


System 
cost C w 


System Average 
Reliability CPU 
R 3 seconds 


Sun et al.[9] 

PFP- G A (this work) 

Big-M-GA(this work) 


2 
7 
7 


(2,1,6,5,0.9396) 
(2,1,6,5,0.9396) 
(2,1,6,5,0.9396) 




0.99992653 9.84 
0.999927 0.2655 
0.999927 0.2590 


Sun et al.[9] 

PFP- G A (this work) 

Big-M-GA(this work) 


4 
9 
9 


(1,1,5,4, .05) 
(2,1,4,4,0.5) 
(2,1,4,4,0.5) 


18.53505 
17.97505 
17.97505 


15.97 

0.3525 

0.2815 



9 Sensitivity Analysis 

To study the performance of our proposed GAs like PFP-GA and Big-M-GA 
based on two different types of penalty techniques, sensitivity analyses (for 
Example-1) have been carried out graphically on the centre of the interval 
valued system reliability with respect to GA parameters like, population size, 
crossover and mutation rates separately keeping the other parameters at their 
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original values. These are shown in Fig.8-Fig.10. From Fig. 8, it is evident 
that in case of PFP-GA, smaller population size gives the better system 
reliability. However, both the GAs are stable when population size exceeds 
the number 30. From Fig. 9, it is observed that the system reliability is stable 
if we consider the crossover rate between the interval (0.65, 0.95) in case of 
PFP-GA. In both GAs, it is stable when crossover rate is greater than 0.8. 
In Fig. 10, sensitivity analyses have been done with respect to mutation rate. 
In both GAs, the value of system reliability be the same. 

10 Conclusions 

In this chapter, the problems of redundancy allocation problems of series 
system, hierarchical series-parallel system, complicated system and reliability 
network system with some resource constraints have been solved. In those sys- 
tems, reliability of each component has been considered as imprecise number 
and this imprecise number has been represented by an interval number which 
is more appropriate representation among other representations like, random 
variable representation with known probability distribution, fuzzy set with 
known fuzzy membership function or fuzzy number. For handling the resource 
constraints, the corresponding problem has been converted to unconstrained 
optimization problem with the help of two different parameter free penalty 
techniques. Therefore, the transformed problem is of unconstrained interval 
valued optimization problem with integer and/or mixed integer variables. 
To solve the transformed problems, two different real coded GA based on 
different fitness functions have been developed for integer and mixed integer 
variables with interval valued fitness function, tournament selection, crossover 
(uniform crossover for integer variables and whole arithmetical crossover for 
floating point variables), mutation (uniform mutation for integer variables 
and boundary mutation for floating point variables) and elitism of size one. 
In the existing penalty function technique, tuning of penalty parameter is a 
formidable task. However, here tuning of parameters is not required as these 
are penalty parameter free techniques. From the performance of GAs, it is 
observed that the GAs with both fitness functions due to different penalty 
techniques take very lesser CPU times with very small generations to solve 
the problems. It is clear from the expression of the system reliability that the 
system reliability is a monotonically increasing function with respect to the 
individual reliabilities of the components. Therefore, there is one optimum 
setup irrespective of the choice of the upper bound or lower bound of the com- 
ponent reliabilities. As a result, the optimum setup obtained from the upper 
bound/lower bound will provide both the upper bound and the lower bound 
of the optimum system reliability. These approaches have wider applicabili- 
ties in solving the constrained optimization problems arisen in every sector 
of real life situation. However, as the proposed techniques are parameter free, 
these do not require the tuning of penalty parameter. 
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Summary. In this chapter, we take advantage of particle swarm optimiza- 
tion to build fuzzy systems automatically for different kinds of problems by 
simply providing the objective function and the problem variables. Parti- 
cle swarm optimization (PSO) is a technique used in complex problems, in- 
cluding multi-objective problems. Fuzzy systems are currently used in many 
kinds of applications, such as control, for their effectiveness and efficiency. 
However, these characteristics depend primarily on the model yield by hu- 
man experts, which may or may not be optimized for the problem. To avoid 
dealing with inconsistent during the fuzzy systems generation, we used some 
known techniques, such as the WM method, to help evolving meaningful 
rules and clustering concepts to generate membership functions. Tests using 
three three-dimensional functions have been carried out and show that the 
evolutionary process is promising. 



1 Introduction 

Fuzzy systems [5J form an important tool to model complex problems based 
on imprecise informations and/or in situations where a precise result is not 
of interest and an approximation is sufficient |16j . The performance of a 
fuzzy system depends on the expert's interpretation, which leads to in the 
generation of the rule base and membership functions of the system. To 
minimize this dependency, some methods are being used in the attempt to 
automatically generate the components required in a fuzzy system. For the 
membership functions, clustering-based algorithms, such as Fuzzy C-means 
and its generalizations as Pre-shaped C-means 0, are usually used. Other 
approaches also exist |10j . The major difficulty in the development of fuzzy 
systems consists of the definition of membership functions and rules that 
provide the desired behavior of these systems. 

Swarm Intelligence is an area of artificial intelligence based on collective 
and decentralized behavior of individuals that interact with each other, as 
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well as with the environment [T] . PSO is a stochastic evolutionary algorithm, 
based on swarm intelligence, that searches for the solution of optimization 
problems in a specific search space and is able to predict the social behavior 
of individuals according to defined objectives [6]. 

Methods based on examples, such as the Wang-Mendel or WM method 
[13] , are usually used for automatic rule generation. Also, there are many 
research works that exploit evolutionary algorithms (EA), both to optimize 
the rule base and the membership functions. In [3], genetic algorithms (GA) 
are used to generate the rule base, with candidate rules pre-selection. In [9], 
the authors use EA to generate fuzzy systems that are more compact and 
more interpretable by humans. In [T3], the authors use clustering techniques 
and GA to define good sets of rules for classification problems. In [5|, the 
authors use evolutionary technique and GA to generate fuzzy systems from 
some given knowledge bases. 

In this paper, we developed an algorithm based on PSO to generate fuzzy 
systems for any kind of problem, provided an objective function that may 
be continuous or discrete. Using simple informations, such as variable names, 
the corresponding domains and the objective function, this algorithm can 
yield an appropriate fuzzy system. Some tests were performed with a known 
control surface to validate the effectiveness of the tool. 

The rest of this paper is organized in five sections. Firstly, in Section [3J 
we explain briefly the principle behind PSO. Then, in Section[2J we describe 
the WM method of rules generation. After that, in Section @J we give details 
about the proposed method for the automatic modeling of fuzzy systems 
using PSO. For this purpose, we first define the structure of a particle and 
the coordinates used to position it within the search space as well as the 
fitness function of the represented system. Then, in Section [SJ we present 
the obtained results to model a commonly used control surface. Last but not 
least, in Sectional we conclude the reported work and give some directions 
for future research. 



2 Particle Swarm Optimization 

During a particle swarm optimization process, each particle is mapped into 
a position in the search space, which n-dimensional. The particle position 
is updated in each iteration. For the position update particle i, the velocity 
related to each of the search space directions is used. The velocity is the 
clement that promotes the movement of the particles and is calculated as in 

© hheiei. 

Vi(t + 1) = wvi(t) + Ciri (£»(*) - Xi(t)) + c 2 r 2 (x(t) - Xi(t)), (1) 

where w is called inertia coefficient, r\ and r 2 are random numbers chosen 
in the interval [0,1], c\ and c 2 are positives constants called as social and 
cognitive coefficients, Xi (t) identifies the best position achieved by the particle 
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in the past and x(t) is the best position, among all the particles, achieved in 
the past. The position of the particle is updated as described in (0). 

Xi (t+l) = Xi(t) + Vi(t+l). (2) 

The velocity guides the optimization process [6] , reflecting both the particle 
experience and the information exchange between particles. The experimental 
knowledge of a particle refers to the cognitive component, which is propor- 
tional to the distance between the particle and its best position, found so far. 
The information exchange between particles refers to the social component 
of the velocity equation (|TJ) . 

To avoid that a particle leaves the search space, it is necessary to use 
a parameter that bounds the velocity [BJ. This is known as the maximum 
velocity v max , it allows a higher granularity on the search control. Therefore, 
before executing the update defined in (J2]), the velocity is analyzed with 
respect to the criterion defined in Q. 

{Vi(t + 1) if Vi(t + 1) < v max 
(3) 
V max if Vi(t + 1) > V max 

In (fTJ), one can observe three terms that interfere in the velocity computa- 
tion 6], which are: 

• The previous velocity, wvi(t), is used to prevent particle i to suffer a drastic 
change in direction. This component also is called of inertia component. 

• The cognitive component, c\r\(xi{£) — Xi(t)), quantifies the performance 
of particle i with respect to previous performances. This component was 
defined by Kennedy and Eberhart as the "nostalgia" of the particle [5] . 

• The social component, C2T2(x(t) — Xi(t)), quantifies the performance of 
particle i with respect to the performance of the set of included particles. 
The effect of this term is to attract the particle to the best position found 
by the particles set. 

The value assigned to each parameter of PSO algorithm is essential in the 
search process evolution. Below are related some set of values considered good. 

• The inertia coefficient, w, controls the relation between exploration and 
exploitation [TJ]. The values near 1 are considered good, but values bigger 
than 1 are not so good so are very small values [Bj. Values bigger than 1 
tend to leave the particles with a very high acceleration, promoting a high 
divergence, while very small values can make the search too slow. 

• The cognitive coefficient, c\, and the social coefficient, C2, yield a better 
performance when these are balanced, i.e., c\ = ci |B]. 

• The factors r\ and r-i define the stochastic nature of the cognitive and 
social contributions. Random values are selected in the range [0,1] [BJ to 
each factors. 
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• The maximum velocity is defined for each of the dimensions of search space 
and can be formulated as a domain percentage [6 ], v max = 5(x rnax ~x m i n ), 
where x max and x m i n are the maximum and minimum domain values 
respectively and 8 is a value in the range [0, 1]. 

• The number of particles define the possibility to cover a range of the search 
space in every iteration of the algorithm. A high number of particles allows 
a better coverage of the search space but requires a considerable computa- 
tional power. Empirical studies show that PSO achieves optimal solutions 
using ten to thirty particles [6]. 

LetpBest be the best position found by a particle and gBest be the best 
position among those found by all the particles. Algorithm Q] is describes 
the PSO optimization process. A given maximal iteration number and the 
predefined fitness values can be used as a stop criterion. 



Algorithm 1. Particle swarm-based optimization algorithm (PSO) 
T 

2 

3 

4 

5 

6 

7 

8 

9 
10 
11 
12 
13 
14 
15 
16 
17 
18 



for i := 1 until total jparticulas do 
Initialize particle i information; 
Initialize random position of particle i; 
Initialize random velocity of particle i; 
end for 
repeat 

for i :— 1 until total jparticulas do 
Calculate fitness of particle i; 
if (fitness better than pBesti) then 

Update pBesti with the new position; 
end if 
if (fitness better than gBest) then 

Update gBest with the new position; 
end if 

Update velocity of particle i; 
Update position of particle i; 
end for 
until (stopcriterion = true) 



The main characteristic of the PSO algorithm is the social interaction 
[5], which makes the individuals able to learn with the group and use the 
acquired knowledge. This allows the particles follow the path that leads to 
more success, as it happens in nature. 

A topology is the way that can be used to understand the interaction 
between a particle and the group in which it is inserted. In practical terms, a 
topology in PSO is defied by the manner the position update in performed. 
Therefore, it is defined by the way the velocity is computed because it through 
the velocity the knowledge about the best path is transmitted between the 
particles of the swarm [5] [TT] . 
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(d) Cluster 



(e) Pyramid (f) Von Neumann 

Fig. 1 Topologies for PSO 



Figure 1 (a) shows the ring topology, also known as Local Best PSO. In this 
topology, each particle is connected to two others, which are in its immediate 



neighborhood. Figure 1(b) shows the star topology, also known as Global Best 



PSO. In this topology, each particle is connected to all those of the swarm, 
i.e., a particle neighborhood is formed by all existing particles. 



Besides these two types of topologies, Figure 1(c) shows all particles con 



nected to central particle that controls the flow of information. In Figure 1(d) 
the particles are clustered and a each particle is considered neighbor of all 
those of its corresponding cluster, yet some of the cluster particle allow for a 



connection with the other in the neighboring clusters. In Figure 1(e), shows 
the pyramid topology which connects particles in triangles. Lastly, Figure 
1(f) | shows the Von Neumann structure, wherein the particles are connected 
via a grid structure [BJ. 



3 Rule Generation Methods 

The rule generation method referred to as Wang-Mendel (WM) [TS] uses 
an input-output data set for the problem at hand, to generate a rule set of 
fuzzy systems. The input-output data set is usually provided as (x p ; y p ), p = 
1,2, ... ,N, wherein x p £ R n and y p £ R. This method extracts the rules 
that best describe how the output variable y £ R is influenced by the n input 
variables x = (x\, ...,x n ) £ R n , based on the provided examples. 
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(c) Output Variable y 
Fig. 2 Example of input-output data set for rules generation 



For instance, assuming two data sets to a system with two input variables 
x\ and X2 and an output y, that are (x\, x\, y 1 ) and (xf, x%, y ), and the 
membership functions showed in the graphics of the Fig. [51 To obtain the 
rules represented by these two sets, first we must get the degree of confidence 
using the membership functions, for each data set. In this case, we have: 

• x\: degree 0.67 in A\ and 0.11 in A2; 

• x\: degree 0.16 in B\ and 0.80 in B2; 

• y 1 : degree 0.66 in CI; 

• x\: degree 0.25 in A2 and 0.68 in A3; 

• x\: degree 0.10 in B3 and 0.58 in B4; 

• y 1 : degree 0.39 in C2 and 0.51 in C3. 

In order to generate the rules, we always keep the membership functions in 
which the variable has the highest degree, and so we discard the functions 
that have lower degree. So, for the data sets defined in Fig. [2], we have (j4|). 
So the first rule would be "If x\ is A\ and xi is B1 then y is CI" and the 
second, "If X\ is A3 and x-i is B4 then y is C3" . 



(x\,xl,y l ) = [si (0.67 in Al), 4(0.80 in B2),y 1 (0.66 in CI)] 
2 ^ - r xf(0.68 in A3),xl(0.58 in B4),y 2 (0.51 in C3)] 



x\^l,y 2 ) 



(4) 
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Note that each data set generates one single rule. Considering a real sys- 
tem, it is very possible that these rules can be conflicting rules. To overcome 
this problem, one can associate degrees of confidence to each generated rule, 
using the degree of relevance of each rule term. Equation [5] shows how this 
degree can be computed: 

C(Rule) = fi(xi) x fi(x 2 ) x fi(y), (5) 

wherein C is degree of Rule and /x(xl), /x(x2) and n(y) are the degree of rele- 
vance of each rule term. In the case of the first rule, the associated confidence 
degree would be C(Rule t ) = 0.67 x 0.80 x 0.66 = 0.35376. 

There are also methods based on genetic algorithms. In [3], the authors 
use the WM method to generate the initial rule set and then apply their own 
genetic algorithm on some classification rate of the rules. This method is only 
used for classification problems. 



4 Proposed Automatic Generation 

In this work, PSO is used to evolve the fuzzy systems parameters of the 
Mamdani type [5], both for rules and membership functions. The search 
algorithm is based on these two elements and always tries to improve the 
solution at hand. However, the functions are not modified in every iteration, 
unlike the rules, whose modification obeys to pre-determined update rate, 
that is defined at the beginning of the evolutionary process. The purpose 
is to maintain the functions stable for some time, giving more time for the 
algorithm to search for more appropriate rules for those functions. At the end 
of each execution, when the algorithm reaches the stop criterion, it returns 
the best solution found. 

There are four important aspect that define the performance of the PSO 
search. These are the particle representation, the position coordinates of a 
given particle in the search space, the fitness function that allows us to de- 
termine how appropriate is the fuzzy system associated with a given particle 
and how to update the system represented by the particle at hand. 

4-1 Representation 

Each particle is associated with a fuzzy system and a position in the search 
space, that is represented by an n-position vector, where n depends on the 
number of used variables. In this work, the fuzzy system is defined by an 
hierarchical structure as described in Fig. [3J 
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Rule Set 

-Variables 

— Membership Functions 



Rules 



Antecedents 



— Consequents 
Fig. 3 Fuzzy representation structure 

4-2 Particle Position and Movement 

The position vector has one entry for the rules number, another one for 
the completeness factor, and m positions for the number of the functions, 
where m is the number of the system variables. Thus, the position dimension 
is dependent on the number of variables of the problem. The completeness 
factor is a criterion that measures the discontinuity between functions in the 
variables domain (see Section I4.3[) . The mutation operator determines how 
the update of each one these items is performed. This update promotes the 
movement of particles on the search space. In this work, we used three kinds 
of mutations: 

• If the velocity relating to the number of rules is positive, then we increase 
the rules number. Otherwise, we decrease it. 

• Changing the number of functions for each variable of the fuzzy system 
follows the same criterion, given above. 

• The change of the completeness is performed increasing the width of a 
function. Thus, the tendency is to reduce the empty space in the domain, if 
any. Similarly, reducing the width of the function, we alter the distinctness 
between the available ones. The more positive the velocity is, the bigger 
the increase in domain of one of the functions. If the more negative the 
velocity is, the smaller the decrease in the domain. 

4-3 Fitness Function 

Inspired by the work reported in [5] , the fitness of each particle is defined as 
in ©: 
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F = -100 x wi x log(MS£) 
+50 x w 2 (l-C r .) 

+50 x w 3 (l - C/) (6) 

+50 x W4 (l - P D ) 
+50xu 5 (l-P c ), 

wherein MSE is the mean-square error of the difference between the returned 
value by the objective function (jjh) and the returned value by the fuzzy model 
(Vh)-> as 0i wherein Njj is the number of data. 

1 N 
MSE=—^2(y h -y^f, (7) 

D h=l 

The C r term represents the relation between the amount of rules presents 
on the model and the total of possible rules, and C/ represents the relation 
between the amount of functions presents on the system and the total of 
possible functions, as in ©: 



c ,■ 



N R 

(8) 



r _ N F 

Term Pjj is a criterion that measures the distinctness between the member- 
ship functions of the variables, defined in ©: 



ft^ENrEo. W 




wherein Ny is the number of variables, N l s is the total possible interval of 
overlap between functions of the i-th variable, A^ is the width of the h-th 
overlap and Xi is the width of the variable domain. In order to the determine 
Xih, it is necessary to define the level £d, drawing a horizontal line, crossing 
all the functions, as showed in the Fig. 4(a)| Term Pc is a criterion that 



represents the completeness of the membership functions, in relation to the 
domain, and is defined as in (fT0|) : 

, N v ( N* D 



wherein N l D is the total possible number of discontinuity between functions 
of the i-th variable, 7^ is the width of h-th discontinuity and Xi is the width 
of the variable domain. 

In order to determine "fih, is necessary to define the level £c> in a similar 
way to £d, as shown Fig. |4(b)[ 
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(a) Overlap 



(b) Discontinuity 



Fig. 4 Overlap and discontinuity illustrations 

The coefficients u>i, u>2, W3, W4 and uj§, control the contribution of each 
term of © in the evaluation of the fuzzy system associated with the particle. 

This evaluation function covers the many required criteria of a fuzzy sys- 
tem. These are the precision, by the error quantification; the compactness, 
by the relation between the number of the rules and functions of the model 
and the total possible number; and the interpretability, by measuring of the 
distinctness and completeness, providing a complete model evaluation. 



5 Results and Tests 

The WM method [15] , introduced in Section [3l was used to initialize rules 
of the fuzzy systems of each particle. Besides, the concept of clustering was 
used in the membership functions generation, to decrease the possibility of 
yielding functions that are incompatible with the fuzzy system. 

Experiments were performed to prove the effectiveness of the technique. 
TableQ]show the used parameters by the PSO-based algorithm. Fig. [5] depicts 
the graphics of the original function z = seno(xy). The best result obtained 



Table 1 Values of the algorithm parameters 



Parameter 


Value 


Social coefficient 


1.5 


Cognitive coefficient 


1.5 


Inertia coefficient 


1 


Number of particles 


20 


Completeness 


0.25 


Number of iteration 


5000 


Total of rules WM 


6 


Kind of function 


Gaussian 


Wi - U>5 


1 
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Fig. 5 Original function z = seno(xy) 




2 -2 



Fig. 6 Best result obtained for function z — seno(xy) 



for the approximation of this function using the fuzzy system generated by 
the PSO-based algorithm is shown in Fig. [5J 

Figure [7] shows the configuration input and output variables of the one of 
the systems yield by the PSO algorithm. The numbering of the membership 
function does not follow the order in which these appear in the domain. This 
due to the Java system used to program the process. 
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Domi'nio da Variavel 
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(b) Variable Y 




(c) Variable Z 

Fig. 7 Membership functions of the fuzzy variables used in the generated system 
for function seno(xy) 



Figure \E\ shows the normalized mean square error (NMSE), that was com- 
puted for each point used to test the validity of the generated fuzzy systems. 
The error is normalized between and 1. 

The general average of the introduced error, considering all the points is 
V = 0.1359. The standard deviation of the error is 8 = 0.0537. 



, 0,15 
' 0,10 




Fig. 8 Error averages for the fuzzy system for function seno(xy) 
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Fig. 9 Original function z = e 



xsin(iry) 




Fig. 10 Best result obtained for function z = e 



xsin(iry) 



Fig. IH1 depicts the graphics of the original function z = e xsln ^ y \ The best 
result obtained for the approximation of this function using the fuzzy system 
generated by the PSO-based algorithm are shown in Fig. [TU1 

Figure ITTI shows the configuration of the input and output variables of the 
resulting fuzzy system. 

Figure [T2] shows the average of the error introduced. In this case, the 
general error average is V = 0.2514, while the standard deviation is 5 = 
0.1364. 
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Dominio da Variavel 

(a) Variable X 




Dominio da Variavel 

(b) Variable Y 




Dominio da Variavel 

(c) Variable Z 

Fig. 11 Membership functions of the fuzzy variables of the system generated for 
the function z = e X3m(ny) 




Fig. 12 Error averages for the fuzzy system for the function z = e 



xsin(iry) 



Table [5] shows a comparison between the results obtained by automatic 
generation of fuzzy systems using genetic algorithms [12] and using the pro- 
posed method using PSO. As before, the error is computed using the Nor- 
malized Mean Square Error (NMSE). 
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Table 2 Comparison ol the performance of GA-based vs. PSO-based fuzzy system 
automatic generation for functions seno(xy) and e xsm (^y) 



Function 


Algorithm 


Fitness 


Deviation of fitness 


Error 


Deviation of error 


seno(xy) 


GA 
PSO 


54 
275.5 


8.0 
3.8 


0.0288 
0.1059 


0.0007 
0.0537 


xsin(ny) 


GA 

PSO 


197 

280.5 


35.0 
16.8 


0.006 
0.176 


0.0032 
0.1656 



6 Conclusion 

In this paper, we illustrated the use of PSO to automatically generate the 
fuzzy rules, fuzzy variables together with the corresponding membership func- 
tions of fuzzy systems. We described the particle representation, its movement 
in the search space and we provided a fitness function that allows us to as- 
sess the appropriateness of the evolved fuzzy system. This experience showed 
that the performance of the evolutionary process is very much dependent on 
the choice of the parameters, such as the number of membership functions 
per variable as well as on the number of rules allowed in the system. More 
tests are being carried out in order to synthesize discrete functions into fuzzy 
systems. 
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Wind energy is an important source of renewable energy, and reliability is a 
critical issue for operating wind energy systems. The Canadian wind energy 
industry has been growing very rapidly. The installed wind energy capacity 
in Canada in 2008 was approximately 2,000 mega watts (MW), which is less 
than one percent of the total electricity. It is believed that wind energy will 
satisfy 20% of Canada's electricity demand by 2025, by adding 55,000MW 
of new generating capacity [1]. Operation and maintenance costs account 
for 25-30% of the wind energy generation cost. Currently, the wind turbine 
manufacturers and operators are gradually changing the maintenance strat- 
egy from time-based preventive maintenance to condition based maintenance 
(CBM) [2-5]. In this article, we review the current research status of main- 
tenance of wind turbine systems, and discuss the applications of artificial 
neural networks (ANN) based health prediction tools in this field. A CBM 
method based on ANN health condition prediction is presented. 



1 Maintenance Optimization of Wind Turbine Systems 

Maintenance management for wind power generation systems aims at re- 
ducing the overall maintenance cost and improving the availability of the 
systems. The existing maintenance methods for wind energy systems can be 
classified into corrective maintenance, preventive maintenance and condition 
based maintenance (CBM) [6]. In corrective maintenance, maintenance ac- 
tivities are performed after failure occurs. There are generally multiple wind 
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turbines in a wind farm, and each wind turbine consists of multiple compo- 
nents, such as gearboxes, main bearings and generators, as shown in Fig. 1 
[7] . Corrective maintenance actions may be performed for a component after 
each failure, or they can be performed after multiple components have failed 
to maintain the multiple failed components simultaneously [1]. Mart?nez et 
al demonstrated there is great need for corrective maintenance optimization 

[5]- 

The preventive maintenance is classified to time-based maintenance and 
age-based maintenance. In time-based preventive maintenance, the mainte- 
nance activities are typically carried out based on the predetermined time 
interval, say every 6 months. Age-based maintenance is employed when the 
component reaches a pre-defined age. However, age-based maintenance is not 
suitable for a multiple wind turbines farm due to expensive fixed cost, which 
is incurred whenever a preventive maintenance is performed once a com- 
ponent reaches the preventive maintenance age values [8]. To further study 
fixed-interval preventive maintenance, Andrawus developed the delay-time 
approach to optimize scheduled inspection plan, and studied a 26?600kw 
wind farm in the case study [9]. 

Condition based maintenance aims at achieving reliable and cost-effective 
operation of engineering systems. In CBM, condition monitoring data, such 
as vibration data, oil analysis data and acoustic data, are collected and pro- 
cessed to determine the equipment health condition; Future health condi- 
tion and thus the remaining useful life (RUL) of the equipment is predicted; 




Fig. 1 Wind turbine components [7] 
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And optimal maintenance actions are scheduled based on the predicted fu- 
ture equipment health condition, so that preventive replacements can be per- 
formed to prevent unexpected failures and minimize total maintenance costs. 
A life cycle cost approach was adopted to evaluate the financial benefit us- 
ing condition monitoring system, a tool for implementing CBM policy [10]. 
In [11] a multi-state Markov decision mechanism was used to estimate the 
wind turbine degradation process based on which the optimal maintenance 
scheme is devised. Tian et al. [8] developed a CBM method for wind turbine 
systems, based on the health condition prediction information obtained from 
ANN prediction models. 

ANN methods have been used to investigate various problems in wind 
turbine systems. Yurdusev et al. proposed a method for determining the op- 
timum tip speed ratio in wind turbines using ANN [12]. Jafarian and Ranjbar 
developed a combined method to estimate annual energy output of a wind 
turbine based on fuzzy modeling techniques and ANN [13]. Intelligent meth- 
ods have also been applied to condition monitoring of wind turbines [14]. 



2 Critical Wind Turbine Components and Their 
Failure Modes 

The critical components of a wind turbine system are discussed in this section. 
Critical wind turbine components include blades, gearboxes, main bearings 
and generators. Wind turbine blades are designed to collect energy from the 
wind and then transmit the rotational energy to the gearbox via the hub 
and main shaft. The number of blades and total area they cover affect wind 
turbine performance. Most wind turbines have only two or three blades on 
their rotors, the reason is that the space between blades should be great 
enough to avoid turbulence, so that one blade will not encounter the dis- 
turbed, weaker air flow caused by the blade which passes before it. In an 
offshore environment, where corrosion is a critical factor to be considered, 
blade material often preferred is corrosion resistant, also the possibility of 
achieving high strength and stiffncss-to-weight ratio. Blades failures include 
cracks arising from fatigue, materials defects accumulating to critical cracks 
and ice build-up are known to cause failure. 

All modern wind turbines have spherical roller bearings as main bearings. 
Main bearing reduces the frictional resistance between the blades, the main 
shaft and the gearbox while it undergoes relative motion. The main bearing 
is mounted in the bearing housing bolted to the main frame. Different types 
of wind turbines vary the quantity of bearings and bearing seats. The main 
bearings ensure that wind turbines withstand high loads during gusts and 
braking. Poor lubrication, wear, pitting, deformation of outer race and rolling 
elements may cause its failures. 

The gearbox is one of the most important and expensive main components 
in the wind turbine. It is placed between main shaft and generator, the task 
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is to increase the slow rotational speed of the rotor blades to the higher 
generator rotation speed. However, the gearbox in the wind turbine does 
not change speed just like a normal car gearbox, it always has the constant 
speed increasing ratio. So if a wind turbine has different operational speeds, 
it is because it has two different sized generators that each one has its own 
different rotation speed, or possibly, one generator has two different stator 
windings. The gearbox is connected to the generator by the coupling. The 
coupling is always a flexible unit made from built-in pieces of rubber, normally 
allowing variations of a few millimeters only. The high speed shaft from the 
gearbox is connected to the generator by means of a coupling. The coupling is 
a flexible unit made from pieces of rubber which allow some slight difference in 
alignment between the gearbox and the generator during normal operation. 
Gearbox failures include poor lubrication, bearings and gear teeth failures 
can cause major failures. 

The generator transforms mechanical energy into electrical energy. The 
blades transfer the kinetic energy from the wind into rotational energy, and 
then generator supply the energy from the wind turbine to the electrical grid. 
Generator produces either alternating current (AC) or direct current (DC), 
and they are available in a large range of output power ratings. The gener- 
ator's rating, or size, is dependent on the length of the wind turbine blades 
because more energy is captured by longer blades. Bearings are the major 
cause of failure of generator. Thus, maintenance is mainly restricted to bear- 
ing lubrication. 

Acronyms 

CBM: condition based maintenance, 
ANN: artificial neural network. 



3 Component Health Condition Prognostics Using 

ANN 

3.1 Health Condition Prediction 

The objective of health condition prognostics is to predict the equipment fu- 
ture health conditions and thus the remaining useful life. At each inspection 
point, the condition monitoring measurements are collected, and the health 
condition prognostics methods can be used to produce the predicted failure 
time value or the remaining useful life value, and some prognostics methods 
can give the associated prediction uncertainties as well. The health condition 
prediction methods can be divided into model-based methods and data-driven 
methods. The model-based methods, also known as the physics-of-failurc 
methods, perform prognostics using equipment physical models and damage 
propagation models. Model-based prognostics methods have been reported 
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for components such as bearings (Marble et al [15]) and gearboxes (Kacprzyn- 
ski et al [16], Li and Lee [17]). The key limitation of the model-based methods 
is that for some components and systems, authentic physics-of-failurc models 
are very difficult to build because equipment damage propagation processes 
and dynamic responses are very complex. Data-driven methods directly uti- 
lize the collected condition monitoring data for health condition prediction, 
and do not require physics-of-failure models. Examples of the data-driven 
methods include the proportional hazards model method developed by Ban- 
jevic et al [18], the Bayesian prognostics methods [19], and the ANN based 
prognostics methods [20-22]. 

Outputs of the prognostics methods are predicted failure time and the 
associated uncertainty. That is, at a certain inspection point, the predicted 
failure time distribution can be obtained for the component being monitored. 
Among various data-driven methods, ANN based methods have been consid- 
ered to be very effective and flexible tools for component health condition 
prognostics. In this work, we use the ANN prediction approach developed in 
[21]. The ANN model used in this approach is shown in Fig. 2, which is a 
feedforward neural network model with one input layer, two hidden layers 
and one output layer. The inputs of the ANN are the component age values 
and the condition monitoring measurements at the current and previous in- 
spection points. In the example of the ANN model shown in Fig. 2, there are 
two condition monitoring measurements. Specifically, U is the age of the com- 
ponent at the current inspection point i, and ti-i is the age at the previous 
inspection point i — 1. z\ and zl_ x are values of measurement 1 at the current 
and previous inspection points, and zfand zf_ 1 are values of measurement 2 
at the current and previous inspection points. The output of the ANN model 
is the life percentage at current inspection time, denoted byP^. For example, 

Input Hidden Output 

Layer Layers Layer 

t 




Fig. 2 Structure of the ANN model for component health condition prediction 

[21] 
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if the failure time of a component is 850 days and the age of the component 
at the current inspection point is 500 days, the life percentage value would 
be P z = 500/850 x 100% = 58.82%. 

The ANN model utilizes suspension histories as well as failure histories. 
A failure history of a unit refers to the period from the beginning of the 
component life to the end of its life, a failure, and the inspection data collected 
during this period. In a suspension history, though, the unit is taken out of 
service before the failure occurs. Once trained, the ANN prediction model can 
be used to predict the remaining life based on the age value of the component 
and the condition monitoring measurements. As mentioned above, the output 
of the ANN model is life percentage, based on which the predicted failure time 
can be calculated. For example, at a certain inspection point, if the age of 
the component is 400 days and the life percentage predicted using ANN is 
80%, the predicted failure time will be 400/80% = 500 days. 

3.2 An ANN Health Condition Prediction Method 

Now we present the ANN health condition prediction method developed by 
Tian ct. al [21]. The procedure of the method is shown in Fig. 3. The detailed 
explanations of the procedure are given in the remainder of this section. 

3.2.1 Constructing the Failure History Training Data Set 

The first step of the approach is to construct the failure history training data 
set, which will be combined with training data set based on the suspension 
histories to train the ANN. Suppose there are J condition monitoring mea- 
surements used in the ANN model. An ANN input vector based on failure 
history / takes the following form: 

1^ = iff, t-l) */, t! Z f, »-l! Z f, t) Z f, t-l! Z f, j) ') Z f, »-l! Z f, i) J (1) 

where tf t i denotes the equipment age at inspection point i in failure history /, 
and z 3 i i represents the measurement j at time tfj. The input vector contains 
the time and the condition monitoring measurements values at the current 
and previous inspection points. The corresponding output value is: 

P t> -#- ( 2 ) 

where TFf represents the failure time for failure history /. Thus, the total 
number of input/output pairs based on the failure histories is: 

F 

N F = ^2(NF f -l), (3) 
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Construct the failure 
history training data set 



Training 



Suspension history 
ID: s = 1 



For suspension history s, specify L discrete 
possible failure time values, and obtain the 
corresponding ANN validation MSE values 



I 



Find the optimal failure time TS S for suspension 
history s based on the ANN validation MSE values 



s = s + 1 



No 




Train the ANN based on the suspension histories 
with optimal failure times and the failure histories 



Prediction 



Predict the life percentage 

using the trained ANN and 

current input values 



Calculate the RUL based on 
the predicted life percentage 




Fig. 3 Procedure of the remaining life prediction approach [21] 
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3.2.2 Finding the Optimal Failure Time for a Suspension History 

The optimal failure time for a suspension history corresponds to the low- 
est validation MSE if we train the ANN using the training set constructed 
based on this suspension history and all the failure histories. For suspension 
history s, we specify L discrete possible failure time values, and obtain the 
corresponding ANN validation MSE values. The discrete failure time values 
are denoted by TSD S ^, TSD S ^ 1 ?, TSD Si l, respectively. These values are 
selected based on the suspension time for the history, TS S . Specifically, we 
can have TSD s j > TS S (1 < I < L) for most of the failure time values, and 
have 1-2 values smaller than TS S , so that we can find the optimal failure 
time based on the validation MSE values at these discrete points. 

For a certain failure time value TSD s j, we can obtain the ANN in- 
put/output pairs for suspension history s. The input vectors take the same 
form as that for failure histories. The ANN output value corresponding to 
the ith. inspection point is given as: 

P s , ^ = rsg^' ( 4 ) 

where t s ^ denotes the equipment age at inspection point i in suspension 
history s. The ANN input/output set includes the input/output pairs based 
on suspension history s and the input/output data set constructed based on 
all the failure histories. Thus, the total number of input/output pairs is: 

F 

N Ss = NS S -1 + J2 ( NF f - 1). ( 5 ) 

/=i 

where NS S represents the total number of inspection points in suspension 
history s. The input/output set is further divided into the ANN training set 
and the ANN validation set: 2/3 of the input/output pairs for the training 
set and 1/3 for the validation set. Specifically, we go through the suspension 
history and the failure histories, and select an input/output pair in every 
three input/output pairs to construct the ANN validation set. The ANN is 
trained using the resilient backpropagation algorithm based on the training 
set and the validation set, and the validation MSE can be obtained. Because 
of the randomness in the training algorithm, typically we cannot obtain the 
exactly same validation MSE value each time. Thus, in this work, we train 
the ANN 30 times, and record the 3 lowest, or best, validation MSE values 
for future use, which are denoted by ve\ r (r = 1, 2, 3). 

The ANN validation MSE values ve s lr (r — 1, 2, 3) can be obtained for 
all the discrete failure time values TSD St i (1 < I < L) for suspension history 
s. Thus, we can obtain totally 3L data points, each containing a validation 
MSE value and the corresponding failure time. In order to find the optimal 
failure time based on the discrete validation MSE values, we need to fit these 
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validation MSE data points. Considering the flexibility and the ability to 
model simple trends, we use the third order polynomial to fit the data points: 

y = ax 3 + bx 2 + ex + d, (6) 

where y represents the validation MSE, x represents failure time, and a, b, c, d 
represent the polynomial coefficients to be determined. The objective of using 
the polynomial function in Equation (J5|) is to build a continuous function 
to represent the change in the validation MSE with respect to the possible 
failure time. Once the polynomial function is obtained, it is easy to find 
the optimal failure time corresponding to the lowest ANN validation MSE, 
using a simple optimization process. The optimal failure time for suspension 
history s is denoted by TS*. To enforce the suspension time constraint, let 
TS* = TS S if TS* is smaller than the suspension time. 

3.2.3 ANN Training Based on the Suspension Histories with 
Optimal Failure Times and the Failure Histories 

Using the procedure in Section 3.3, we can find the optimal failure times 
for all the suspension histories: TS^ (1 < s < S). Now we can train the 
ANN for remaining useful prediction based on the suspension histories with 
optimal failure times and the failure histories. The ANN output value for an 
input/output pair from a failure history is given by Equation (U]), and that 
from a suspension history is given as follows: 

Thus, the total number of input/output pairs is: 

S F 

JVio = E ( NS * ~ 1) + E ( NF f ~ ^' ( 8 ) 

8=1 / = 1 

The ANN training set includes 2/3 of the input/output pairs, and the ANN 
validation set includes the remaining 1/3 of the input/output pairs. Similarly, 
we train the ANN 30 times using the resilient backpropagation algorithm, and 
save the ANN with the smallest validation MSE. 



3.2.4 Remaining Life Prediction Using Trained ANN 

Once the ANN is trained, as discussed in the previous section, it can be 
used for RUL prediction for other equipments being monitored. The age 
and condition monitoring measurements at the current and previous data 
points are used as inputs to the trained ANN, and the current life percentage 
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can be obtained. The RUL is obtained by dividing the current age by the 
predicted life percentage. When new condition monitoring data is available, 
the prediction will be performed again and the RUL will be updated. The 
remaining useful life prediction process stops when the equipment fails or 
when it is preventively taken out of service. 

3.3 Quantification of the ANN Health Condition 
Prediction Uncertainty 

In this section, we present a method for estimating the predicted failure time 
distribution based on the ANN life percentage prediction errors obtained 
during the ANN training and testing processes [23]. 

In the ANN training process, the ANN model is trained based on the 
available failure histories and suspension histories. The ANN model inputs 
include the age data and the condition monitoring measurements at the cur- 
rent and previous inspection points. The output of the ANN model is the 
life percentage of the inspected component at the current inspection point, 
denoted by P,. In the training process, the weights and the bias values of 
the ANN model are adjusted to minimize the error between the ANN output 
and the actual life percentage. After ANN training is completed, the predic- 
tion performance of the trained ANN model is tested using testing histories 
which are not used in the training process. Here, the ANN prediction error 
is defined as the difference between the ANN life percentage output and the 
actual life percentage value at an inspection point in the test histories. That 
is, the ANN prediction error at inspection point k in a test history is equal to 
(Pfc-Pfc), where P& denotes the actual life percentage and P& is the predicted 
life percentage using ANN. Since a test history contains many inspection 
points, with several test histories, we can obtain a set of ANN life percentage 
prediction error values. 

In this study, it is assumed that the ANN life percentage prediction error 
is normally distributed, since the prediction uncertainty is mainly due to 
the capability of the ANN prediction model. With the obtained set of ANN 
prediction error values, we can estimate the mean /i p and standard deviationcr p 
of the ANN life percentage prediction error. Suppose at a certain inspection 
point, the ANN life percentage output is Pt, then the mean of the predicted 
life percentage should be P t — n P , and the standard deviation is still a p . If 
the age of the component at the current inspection point is t, the predicted 
failure time will be tj (P t — n P ), and the standard deviation of the predicted 
failure time will be a p ■ tj (Pt — /j, p ). That is, the predicted failure time T p 
follows the following normal distribution: 

T p ~ N (tj (P t - n p ), a p ■ t/ (P t - n p )) . (9) 
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4 A CBM Approach for Wind Power Generation 
Systems 

In this section, the CBM method developed by Tian et al. [8] for wind power 
generation systems is presented. Suppose there are N wind turbines in the 
wind farm, and we consider M critical components for each turbine. In this 
work, it is assumed that all the wind turbines under consideration are iden- 
tical, and the degradation processes of the wind turbine components are 
mutually independent. 

4-1 Failure Probability Estimation at the Component 
and Turbine Levels 

At the wind turbine component level, condition monitoring data, such as vi- 
bration data and acoustic emission data, can be collected, and failure time 
distribution can be predicted for each component using the prognostics meth- 
ods presented in Section 2. It is assumed that the predicted failure time 
follows the normal distribution, as discussed in Section 2. The failure proba- 
bilities for the wind turbine components, which will be defined later, can be 
calculated based on the predicted failure time distributions, and the CBM de- 
cisions will be made based on the failure probabilities. The failure probability 
for component m in turbine n is defined as follows [23]: 

Pr = ~ "^ lfx _ t , 2 (10) 

J t a\/2TT 

where t is the age of the component at the current inspection point, t p is 
the predicted failure time using ANN, and a is the standard deviation of the 
predicted failure time distribution. Based on the discussions in Section 2, we 
have the following relationships: 

t p = t/ (P t - fX p ),<J = cr p -t/ (P t -Up). (11) 



L in Equation (jlOD is the maintenance lead time, which is defined as the 
interval between the time maintenance decision is made and the time when 
the maintenance actions are performed. The lead time consists of the time 
required to assemble the maintenance team, order the replacement parts, 
prepare the maintenance equipments to perform the maintenance, and travel 
to the wind farm, etc. Thus, the maintenance decisions made at the current 
inspection point can affect the wind turbines only when the lead time has 
passed, and we have no influence on the failures during the lead time. So, it 
is reasonable to base the maintenance decisions on the failure probabilities 
during the maintenance lead time L in order to reduce the failure risks. To 
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reasonably simplify the problem, we assume L is the same for all maintenance 
actions in this studyto simplify our discussion. 

If we focus on the critical components in a turbine, such as rotor, main 
bearing, gearbox, generator, etc., the turbine can be treated as a scries sys- 
tem. That is, the failure of any turbine component will cause the failure of 
the turbine. Thus, the failure probability of wind turbine n during lead time 
L can be calculated as follows: 

M 



Pr = 1 - TT ( 1 - Pr ) (12) 

n. A A \ n rn. I 



4.2 The Proposed CBM Policy 

For the purpose of simplifying the descriptions, we use replacement to refer 
to a maintenance action, such as the replacement of the main bearing, or the 
replacement of a faulty gear within the gearbox. Suppose wind turbine com- 
ponents are continuously monitored. Maintenance decisions are made based 
on the failure probabilities of the components and the wind turbines, which 
can be calculated based on the component health condition monitoring and 
prognostics information. 

The proposed CBM policy for the wind power generation systems is sum- 
marized as follows: 

1. Perform failure replacement if a component fails. The maintenance equip- 
ments and replacement parts will be scheduled, and the maintenance team 
will be sent to the wind farm. 

2. Send a maintenance team to the wind farm and perform preventive replace- 
ments if any wind turbine in the wind farm is determined to be maintained. 

3. Perform preventive replacements on components in wind turbine n if 
Pr n > d\, where Pr n is the failure probability of the wind turbine n, 
and d\ is the pre-specificd level 1 failure probability threshold value at the 
turbine level. 

4. If turbine n is to be performed preventive replacement on, perform pre- 
ventive replacement on its components in order to bring the turbine failure 
probability down to below efe. d-i is called the level 2 failure probability 
threshold value at the turbine level. 

As can be seen, once the two failure probability threshold values, d\ and di-, 
are specified, the CBM policy is determined. 

4-3 CBM Optimization Model and Solution Method 

Based on the above CBM policy, the CBM optimization model can be briefly 
formulated as follows: 
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min Ce (di, d 2 ) 

s.t. (13) 

< d 2 < di < 1 

where Ce is the total expected maintenance cost per unit of time for a certain 
CBM policy defined by the two failure probability threshold values d\ and 
d 2 . d\ and d 2 take real values between and 1, and d 2 < d\. The objective 
of the CBM optimization to find the optimal d\ and d 2 values to minimize 
the total maintenance cost. The optimization functions built in Matlab can 
be used to solve this optimization problem, and find the optimal threshold 
failure probability values. 

Before the optimization can be performed, we need to first be able to 
calculate the cost value Ce given two failure probability threshold values d\ 
and d 2 . Due to the complexity of the problem, it is very difficult to develop a 
numerical algorithm for the cost evaluation of the CBM policy for the wind 
power generation systems. In this paper, we present a simulation method 
for the cost evaluation. The flow chart for the procedure of the simulation 
method is presented in Fig. 4, and detailed explanations of the procedure are 
given in the following paragraphs. 

Step 1: Building the ANN prediction model. For each type of wind turbine 
component, determine the life time distribution based on the available failure 
and suspension data. Weibull distributions are assumed to be appropriate 
for components lifetime, and the distribution parameters a m and j3 m can be 
estimated for each component m. For each type of component, based on the 
available failure and suspension histories, an ANN prediction model can be 
trained, and the mean and standard deviation of the ANN life percentage 
prediction error, which are denoted by fJ, PtTn and a PtTn , respectively, can be 
calculated. 

Step 2: Simulation initialization. As mentioned earlier, suppose there are 
N wind turbines considered in the wind farm, and M critical components 
are considered for each turbine. Specify the maximum simulation time Tjyiax, 
and the inspection interval Tj. Tj can be set to be a small value, say 1 day, 
so that we can approximately achieve continuous monitoring. Or we can set 
Tj to be a bigger value, say 10 days, to improve computation efficiency and 
achieve reasonably accurate results. For each component to, specify the cost 
values, including the failure replacement cost c/ iTO and the variable preventive 
replacement cost c P]Tra . The fixed cost of maintaining a certain wind turbine, 
c p ,t, and the fixed cost of sending a maintenance team to the wind farm, 
CFarm, also need to be specified. The total replacement cost is set to be 
Ct = 0, and will be updated during the simulation process. At time set 
tABS = 0j generate the real failure times for each component in each turbine. 
That is, for component to in turbine n, generate a real failure time TL n _ m 
by sampling the Weibull distribution for component m with parameters a m 
and (3 m . Thus, at time 0, the age values for all the components are 0. 
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Step 1: Building the ANN 
prediction model 



Step 2: Simulation initialization 



Step 3 : Component health 

condition prognostics and 

failure probability calculation 



Step 4: CBM decision making, 

cost update, and component age 

and real failure time values update 




Step 5 : Calculating the total 
expected replacement cost Ce 



Fig. 4 Flow chart for the proposed simulation method for cost evaluation 

Step 3: Component health condition prognostics and failure probability cal- 
culation. At a certain inspection point when the time is £abs>0, the age 
of component m in turbine n is represented by i n , m , and its real failure 
time is known at this point, which is TL nrn . For each component, gener- 
ate the predicted failure time, TP n ^ m , by sampling the normal distribution 
N (TL nym , <j p ■ TL nym ). Based on the discussion in Section 2, the predicted 
failure time distribution can be obtained as N (TP„ im , a p -TP nym ). Thus, 
the current failure probability during the lead time for the component is: 



Pr = 

n.rn 



f t„„ 

Jt„ r. 



i+L 



<j p TP n , 



/2-it 



dx 



SI 



(14) 



<r„TP„ 



/2tt 



dx 



The failure probabilities of each turbine can be calculated using Equation 
(|12p based on the failure probabilities of its components. 

Step 4: CBM decision making and cost update. At the current inspection 
point tABS> the CBM decisions can be made according to the CBM policy, 
described in Section 3.2, based on the failure probabilities of the turbines and 
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their components. If the current time tABS has not exceeded the maximum 
simulation time Tm&x, repeat Step 3 and Step 4. 

Step 5: Total replacement cost calculation. When the maximum simulation 
time is reached, that is, tABS = 2\iax, the simulation process is completed. 
The total replacement cost for the wind farm can be calculated as: 



Ce — CV/Tiyiax- 

And the total replacement cost for each turbine is: 



Cet = 



N ■ Tmelx 



(15) 



(16) 



5 An Example 

In this section, an example is used to demonstrate the proposed CBM ap- 
proach for wind power generation systems [8]. Consider a group of 5 wind 
turbines, produced and maintained by a certain manufacturer, in a wind farm 
at a remote site. To simplify our discussion, in this example, we study 4 key 
components in each wind turbine: the rotor (including the blades), the main 
bearing, the gearbox and the generator, as shown in Fig. 5 [24]. 

Assume the Weibull distributions are appropriate to describe the compo- 
nent failure times, and the Weibull parameters are given in Table 1. The 
component lifetime distribution parameters are specified based on the data 
given in Ref. [25] and [26]. The cost data are given in Table 2, including 
the failure replacement costs for the components, the fixed and variable 



Rotor 




Generator 



Gearbox 
Main Bearing 
Fig. 5 Key wind turbine components considered in the example [24] 
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Table 1 Wcibull failure time distribution parameters for major components 



Component 


Scale parameter 

a(days) 


Shape parameter 

13 


Rotor 


3,000 


3.0 


Main bearing 


3,750 


2.0 


Gearbox 


2,400 


3.0 


Generator 


3,300 


2 



preventive replacement costs and the cost of sending a maintenance team 
to the wind farm. The cost data are specified based on the cost related data 
given in Rcf. [1] and [27]. The ANN prediction method can be used to predict 
the failure time distributions of the wind turbine components, and suppose 
the standard deviations of the ANN life percentage prediction errors are 0.12, 
0.10, 0.10, and 0.12, respectively, as shown in Table 3. The standard deviation 
values are selected by referring to that estimated using the bearing degrada- 
tion data in Ref [18] and [28]. The maintenance lead time is assumed to be 
30 days, and the inspection interval is set at 10 days. 

Table 2 Failure replacement and preventive maintenance costs for major 
components 



Component 


Failure 
replacement 
cost ($1000) 


Variable 
preventive 
maintenance 
cost ($1000) 


Fixed 
preventive 
maintenance 
cost ($1000) 


Fixed cost to 
the wind farm 
($1000) 


Rotor 


112 


28 




50 


Main 
bearing 


60 


15 


25 




Gearbox 


152 


38 






Generator 


100 


25 







The total maintenance cost can be evaluated using the proposed simulation 
method presented in Section 3.3. The cost versus failure probability threshold 
values plot is given in Fig. 6, where the failure probability threshold values 
are given in the logarithm scale. It is found that the total maintenance cost is 
affected by the two failure probability threshold values, and the optimal CBM 
policy exists which corresponds to the lowest cost. Optimization is performed, 
and the optimal CBM policy with respect to the lowest total maintenance cost 
can be obtained. The obtained optimal threshold failure probability values 
are: d\ — 0.1585, e?2 = 3.4145?10~ 6 , and the optimal expected maintenance 
cost per unit of time is 577.08 $/day. 
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Table 3 ANN life percentage prediction error standard deviation values for major 
components 



Component 


Standard deviation 


Rotor 


0.12 


Main bearing 


0.10 


Gearbox 


0.12 


Generator 


0.10 



1200 
1000 

03 

§ 800 



O 600 



400 





-30 -4 



L °9(^) ~~ ' Logfd,) 

Fig. 6 Cost versus failure probability threshold values in the logarithm scale 



6 Conclusions 

Wind energy is an important source of renewable energy, and reliability is a 
critical issue for operating wind energy systems. Maintenance optimization 
approaches aiming at improving wind turbine system reliability and reducing 
the overall operating costs. Currently, the wind turbine manufacturers and 
operators are gradually changing the maintenance strategy from time-based 
preventive maintenance to condition based maintenance. Given the failure 
probabilities for components and the system, optimal CBM decisions can 
be made on target wind turbines to be maintained, the maintenance sched- 
ule, and key components to be inspected and fixed. Future research includes 
developing wind turbine component health monitoring methods by utiliz- 
ing data collected under different conditions, building more accurate wind 
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turbine system reliability models considering the degraded system perfor- 
mance, and maintenance optimization for farms with heterogeneous wind 
turbines, etc. 
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Economic dispatch is an important problem in power systems. This chapter 
presents how a method of stochastic optimization, a metaheuristic known as 
CLONALG (CLONal selection ALGorithm), can be applied to the economic 
dispatch problem. The objective function used in the optimization is based 
on Karush-Kuhn- Tucker conditions, thus, guaranteeing a convergence to the 
global optimum. Examples and results are presented showing the method is 
capable of finding the optimum solution while respecting power generation 
limits. 



1 Fundamentals of CLONal Selection ALGorithm 
(CLONALG) 

Clonal selection theory [5] is an widely accepted model for acquired immu- 
nity. According to the theory, the continuous encounter of lymphocytes and 
antigens causes biological triggers that lead to a selection process of the 
lymphocytes, causing the antibody population to have a greater affinity for 
that antigen. This process is, effectively, a minimization of the error on the 
function of antibody encounters with antigens. This theory gave origin, in 
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computing, to the class of clonal selection algorithms of immune artificial 
systems (IAS). 

The theory of clonal selection has the following characteristics: mainte- 
nance of a memory set, selection and cloning of the most activated antibodies, 
death of non-activated antibodies, affinity maturation through the mutation 
process, reselection of clones proportional to the antigen affinity, generation 
and maintenance of diversity. 

Aiming to create a minimization or classification algorithm with similar 
characteristics, the clonal selection algorithm (CSA) was proposed in [6] and 
later, with a more developed concept, the new concept of the algorithm was 
published by [7J and named clonal selection algorithm (CLONALG). The 
pseudo code of the algorithm is detailed in Table [1] 

Table 1 Pseudo code of CLONALG 

[f* , /, A b ] = CLONALG(F obj , N abs , N ael , N new ,p,p, N gens ) 

[At] = generate(Nabs); {Generate initial population} 
for i — 1 : N gens do 

[/] = evaluate(F bj,Ab); {Determine population affinity} 

[Ab„ , f n ] — select(A b , f, N 3e i); {Select best} 

[C,f c ] = clone(Ab n , f n , P); {Proportional cloning} 

[Cmut] — hyper mutate(C, f c , p); {Affinity maturation} 

[/mut] = evaluate(F bj,C m ut); {Determine clone affinity} 
[Ab n ,fn] — select(C m ut, fmut, N sei ); {Select best} 

[A;,,/] = insert(Ab, Ab n , f, fn)', {Insert clones} 

[Ab d ] = generate(N new ); {Generate new antibodies} 

[Ab] = replace(Ab, Ab d , f); {Replace worst} 
end for 

[/] = evaluate(F bj , Ab); {Determine final solutions} 

[/*] = min(f); {Determine best solution} 

The algorithm in it's original version presents affinity proportional cloning, 
similar to the biological inspiration. This is achieved by sorting the antibody 
population by affinity to the antigen in ascending order and then calculating 
the number of clones for each antibody according to: 

N c i nes = round f r- 2 -^- + 0.5 j , (1) 

where i is the antibody rank, i e [l,N abs ]. However, in the optimization 
version of CLONALG the affinity proportional cloning is not effective [TO] , 
The recommended approach is to use the equation given by [7] : 

N c i nes = round {(3 ■ N abs ) . (2) 

The mutation, also known as affinity maturation, is given with a rate 
inversely proportional to the fitness of each antibody. That is, antibodies 
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with a lower affinity to the antigen will suffer more mutations while the 
antibodies with a better affinity will suffer less mutations. This is done by 
first normalizing the objective function of each antibody, /j, in the range 
[0, 1]. Then a random perturbation vector is generated for each clone created 
from the antibodies, this vector is scaled by the mutation rate, which is given 
by the following equation: 

a = eS-P-M . (3) 



1 . 1 Parameters 

CLONALG 's performance depends heavily on the user defined parameters, 
which are: antibody population size, selection pool size, remainder replace- 
ment pool size, clonal factor and mutation factor. Basic output parameters 
are: F bj , the objective function to be minimized; /*, the best objective 
function value found; /, a vector of final values for the objective function; 
and Ab, a matrix of final solutions. 

Antibody population size (iV a b s ) — The total amount of antibodies to 
be used by the algorithm. Antibodies are an analogy to solution vectors for 
an objective function. 

Selection pool size (N se i) — The number of best antibodies to be se- 
lected by the algorithm for the cloning stage. This is similar to an elitism 
approach, like the one used in Genetic Algorithms (GA). Only clones and 
new antibodies are passed on to newer generations, thus, the value of N se i 
determines the selective pressure on the population. Greater values soften 
the selective pressure, ensuring more members of the population are cloned 
and the search is broadened across the search space, causing the population 
diversity to increase. On the other hand, lower values can increase the se- 
lective pressure by allowing only a few solutions with greater affinity to be 
cloned, causing a portion of the antibody population to be dominated by 
these solutions, reducing diversity. The trade-off between exploitation, local 
search, and exploration, global search, is mainly done by the selection of this 
parameter. 

Remainder replacement pool size (N new ) — The number of lower 
affinity antibodies to be replaced by random antibodies on each generation. 
Essentially, this is a mechanism to ensure a minimum, as in additional, diver- 
sity is maintained in the antibody population. In practice, this mechanism 
is only useful when the algorithm is stuck on local optima. This additional 
diversity can be disabled by setting the value to zero. 

Clonal factor (/3) — A scaling factor for the number of clones created for 
each selected antibody. Common values are j3 £ (0,1]. Lower values decrease 
the clone population while higher values increase it. 

Affinity maturation factor (p) — A scaling factor for the affinity matu- 
ration. Since the affinity maturation is inversely proportional to the affinity of 
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antibodies, so is the scaling factor. Higher values decrease mutation diversity 
while smaller values increase mutation diversity. 

2 Economic Dispatch of Electrical Energy 

The economic dispatch problem (EDP) is an important power systems op- 
timization problem with the goal of attaining optimum power dispatch for 
generators while respecting certain restrictions [31] . Conventional methods 
for the dispatch problem use Lagrangian multipliers, which take the prob- 
lem's derivatives into account, yet, the generating unit's characteristics are 
inherently non-linear, creating multiple local minima in the cost function |27| , 
[2"fi] . |34) . In fact, the problem is complex and highly non-linear [5], it is highly 
multi- modal in the prohibited zones and the valve-point loading versions |23| . 
The economic dispatch problem, it's variants and different solution methods 
have been extensively studied with several works like [19], [11], [30], [28], [29] . 
Ql, 0, ESI, 0. B2, US, EB, 0, tonameafew. 

Economic dispatch is the problem of deciding the most efficient, low-cost, 
configuration of a power system, while operational constraints, by dispatch- 
ing the available power generation to meet a certain system load. It's de- 
scribed as a minimization problem of the total cost of power generation in 
a system while satisfying constraints, power limits, on the system's genera- 
tion resources. There are two basic classes of constraints in the problem: i) 
inequality restrictions, where each generator has it's power restricted to min- 
imum and maximum limits; ii) equality restriction, where the total generated 
power must be equal to the system load, or in the version with losses, must 
be equal to the system load plus losses. In mathematical terms, the simple 
version of the problem, which ignores prohibited zones and/or valve-point 
loading effects, is described as: 



Ng 

min ^Ci(P gi ) (4) 

Ng 

subject to: 2_, ^9i = Pd 

4=1 

Min(Pg t ) < Pg t < Max{Pg t ), 

where Ng is the number of generators, i G [1, Ng], Ci is the cost equation 
($/h) of each generator, Pg t is the active generated power (MW) of each 
generator and Pd is the active power demand of the system. The cost equation 
Ci may assume several forms, however, it's most common in the literature in 
the quadratic matrix equation formxomo: 

Ci{x) — x ■ Q ■ x + c ■ x + a, (5) 
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where x is a vector with active generated power of each generator, Q is a 
positive-definite symmetric matrix which defines the second order elements, 
c is a column vector which defines the first order elements and a a column 
vector which defines the constants. Being a minimization problem, methods 
based on derivatives of the cost equations will ignore the a term. 

Being a non-linear problem, iterative methods are necessary for the so- 
lution. Common iterative methods in the literature include the lambda 
method [32]; gradient projection [13]; linear programming [55j; dynamic pro- 
gramming [25]; primal-dual interior points method [22], [14] and predictor- 
corrector method [33 . However, it's a characteristic of non-linear problems 
for the computational cost to increase expressively according to the problem 
dimension [5]. 

The rest of this work is organized as follows. Section [3] describes the eco- 
nomic dispatch problem as an optimization problem, detailing the objective 
function and search space. Section [4] contains examples and presents the re- 
sults of using CLONALG to the minimization of the KKT point error. Finally, 
the conclusion is given in section [5] 

3 Optimization Problem 

Since the economic dispatch problem is essentially defined by Equation^ the 
search space is easily defined as a vector of generated power of each generator, 
that is: 

search space = [Pgi, ...,PgNg}- (6) 

Many metaheuristics, such as CLONALG, can handle lower and upper lim- 
its in the search space. These are defined in the economic dispatch problem 
as the inequality constraints, thus, it's possible to take advantage of CLON- 
ALG's handling of search space limits by removing the inequality constraints 
from the problem, changing the problem definition from Equation 0] to: 



mm 



Ng 

fin J^CiiPgi) (7) 

Ng 

subject to: 2_\ Pg% = Pd- 



i=l 



Removing the inequality constraints from the problem definition decreases 
the problem complexity, leaving only the equality constraint. A common ap- 
proach for metaheuristics to the economic dispatch problem is to minimize 
the total generator costs while considering the power demand constraint with 
a high valued penalty factor, such as: 
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Ng /Ng \ 

F obj =J2 C * ( p 9i) + penalty • ^ P 9l = Pd . (8) 

The problem with this approach is it lacks a guarantee of a global mini- 
mum. Even if the optimization algorithm, whichever is used, stalls indefinitely 
at a minimum, there is no numeric guarantee the global minimum has been 
reached. A more clever approach, which guarantees convergence to the global 
minimum or at least a measurement of how far away from it, is to seek the 
Karush-Kuhn-Tucker (KKT) point Q2|, [18], A x C(x*,X*) = 0, as it has been 
done in ISO]. 



3.1 KKT Point Minimization 

Recalling Equation [7] and expanding it, the following is obtained: 

Ng 

mm J]PgT.Q i .Pg i + cf-Pg i (9) 

i=\ 

Ng 

subject to: %. P& = Pd, 

then utilizing the Lagrangian method to reach the minima: 



Ng / Ng 

C[p gi , X) = Y, ( p 9l ■ ft • Pm + cf • Pgi) +x T [Pd-Y, Pgi 

4=1 \ i=\ 



(10) 



The global minimum is the set of generated power of each generator, Pgi, 
and the incremental cost, A, which satisfies the partial derivatives of Equation 
ITTJ1 in the origin. The new search space is then defined by: 

search space — [Pgi,...,PgN g ,\]. (11) 

Finally, the objective function for the problem is the error to the KKT 
point, given a set of power generation units and an incremental cost: 

_ v-~* f dC \ dC .„„. 

F * = g(^) + ^' (12) 

For a system with a single generator, the objective function can be visual- 
ized in Figure [U where it's easy to see how the global minimum is at the end 
of a long, narrow, almost flat valley. Finding the minimum region, the valley, 
is trivial, but to converge to the global minimum is, however, very difficult. 
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Fig. 1 Objective function for a system with a single active power generator, in 
logarithmic scale 



4 Examples 

The examples presented in this section do not take into account prohibited 
zones or valve-point loading effects. The maturation process of CLONALG 
needs to be slightly modified to increase performance in the economic dispatch 
problem. Considering the narrow valley of the objective function, a Gaussian 
distribution is more useful in the mutation process to increase local search. 
Also, to ensure potential solutions aren't lost, maturation only takes place in 
50% of dimensions in all clones. 

Results are presented for 50 runs of CLONALG for each example with 
the following parameters: antibody population size, N a }, s = 10; selection pool 
size, N se i = 10; remainder replacement pool size, N new = 0; clonal factor, 
(3=1; affinity maturation factor, p = 0.25; maximum generation, N gens ) = 
50000. 

To assess the efficiency of CLONALG in the minimization to the KKT 
point error, two case studies, with and without power loss, of economic dis- 
patch problems are analyzed. 



4-1 Case Study 1 

This case study, presented in Table [3J consists of three generator units with 
configuration given in Tabled 

The convergence results in Tableware presented for running CLONALG 
with 50000 generations, for 50 runs, showing CLONALG's capability to find 
the global optimum. Results for the objective function value, the error to 
the KKT point, and for the generation cost are also presented, showing min- 
imum, mean, standard deviation and maximum values. Table 0] shows the 
best solution found by CLONALG for this case study. 



NO 
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Table 2 Data of case 1 of economic dispatch problem 



G 


a 


b 


c 


Pg min (MW) 


Pg ma , (MW) 


Pd (MW) 


1 


0.1562 


7.92 


561 


150 


600 


850 


2 


0.00194 


7.85 


310 


100 


400 


3 


0.00482 


7.97 


78 


50 


200 



Table 3 Minimization results for the case 1 



Min. F obj 


Mean F tj 


Std. Dev. of Fobj 


Max. F t,j 


5.4502e-08 


6.7467c-06 


1.3786e-05 


6.6261e-05 


Min. Cost ($/h) 


Mean Cost ($/h) 


Std. Dev. Cost ($/h) 


Max. Cost ($/h) 


8194.35610840853 


8194.35612105488 


2.0573e-006 


8194.35612491201 



Table 4 Best solution obtained for the case 1 in 50 runs 



± obj 


Pgt (MW) 


Pg* 2 (MW)jPp 3 * (MW) 


Incremental cost A* ($/MWh) 


5.4502e-08 


393.169843 


334.603752J 122.226405 


9.148263 



4.2 Case Study 2 

This case study, presented in Table [S] is similar to Case 1 , the difference lies 
in the power loss dependent on generator 1. The power loss changes the KKT 
conditions by affecting the Lagrangian derivative of the first generator and of 
the incremental cost, presenting a new challenge to CLONALG's convergence. 

Table 5 Data of case 2 of economic dispatch problem 



G 


a 


b 


c 


Pgrmn (MW) 


Pg max (MW) 


Pd (MW) 


1 


0.1562 


7.92 


561 


150 


600 


850 + 0.01359-Ppi 


2 


0.00194 


7.85 


310 


100 


400 


3 


0.00482 


7.97 


78 


50 


200 



The convergence results in Table [6] are presented for running CLONALG 
with 50000 generations, for 50 runs, showing CLONALG's capability to find 
the global optimum in a different scenario, where a power loss element based 
on generator 1 is introduced in the problem. Results for the objective function 
value, the error to the KKT point, and for the generation cost are also pre- 
sented, showing minimum, mean, standard deviation and maximum values. 
Table [7] shows the best solution found by CLONALG for this case study. 
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Table 6 Minimization results for the case 2 



Min. F obj 


Mean F a bj 


Std. Dev. of F ob:j 


Max. F bj 


3.7424e-08 


0.0012583 


0.0085057 


0.060162 


Min. Cost ($/h) 


Mean Cost ($/h) 


Std. Dev. Cost ($/h) 


Max. Cost ($/h) 


8242.24097636649 


8242.24667247872 


0.0401672683253928 


8242.52501712881 



Table 7 Best solution obtained for the case 2 in 50 runs 



obi 



Pgt (MW) Pg* 2 (MW) 



Pgl (MW) Incremental cost A* ($/MWh) 



3.7424e-08 



374.299972 351.685238 



129.101528 



9.214539 



5 Conclusion 

This work analyzed the usage of CLONALG on the economic dispatch prob- 
lem of electrical energy. The optimization problem was formulated as the min- 
imization of the KKT point error, considering the issue that most problems in 
stochastic minimization have no guarantee of a global optimum. This formu- 
lation may not assure a convergence to the global optimum with stochastic 
optimization algorithms but, at least, provides the knowledge of how far away 
the solution is from the optimum and the objective function is shown to pos- 
sess an easy to find minima region, although the global minimum is hard to 
find due to the valley created by the incremental cost. Minor modifications 
have been made to improve CLONALG 's performance on the problem, such 
as changing the perturbation vector, the mutated clones, from an uniform to 
a Gaussian distribution with mean and variance 1. In addition, a selection 
of dimensions similar to the binary selection used in differential evolution was 
used, where only 50% of the dimensions are effectivelly mutated, or matu- 
rated, reducing losses of potential solutions. The approach was tested on two 
case studies, where case study 1 is a simple power system with three genera- 
tors and case study 2 is an extension of the first, with the addition of power 
losses, to test the robustness of the method. Empirically, results demonstrate 
the approach is successful, showing CLONALG is able to reach the global 
optimum on both case studies. 
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Several approaches for solving multi-objective optimization problems entail a 
form of scalarization of the objectives. This chapter proposes a study of dif- 
ferent dynamic objectives aggregation methods in the context of evolutionary 
algorithms. These methods are mainly based on both weighted sum aggre- 
gations and curvature variations. Since the incorporation of chaotic rules or 
behaviour in population-based optimization algorithms has been shown to 
possibly enhance their searching ability, this study proposes to introduce and 
evaluate also some chaotic rules in the dynamic weights generation process. A 
comparison analysis is presented on the basis of a campaign of computational 
experiments on a set of benchmark problems from the literature. 



1 Introduction 

A multi-objective or vector optimization problem with m > 2 objectives or 
criteria can be stated as follows: 

mmf(x) = (h(x),...,f m (x)y , (1) 

where fj : R" — > M, for j = 1, .. . ,m; X C R" is the set of the feasible 
decision vectors. Generally, the feasible set is implicitly given through a set 
of constraints. Z :— f(X) is the set of all values assumed by the objective 
function on the feasible set; it is a subset of the objective space R m and the 
vectors belonging to Z are called feasible criteria vectors. 
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While a multi-objective problem involving independent design functions 
can be solved by simply minimizing m scalar objective functions separately, 
when the objectives are competitive it is very difficult to obtain a single 
decision vector which minimizes all criteria simultaneously. 

Generally, the solution of a multi-objective optimization problem is a set of 
vectors, the Pareto solutions, that provide a trade-off among criteria. The goal 
of a multi-objective optimizer is to achieve a set of Pareto optimal solutions. 
Since every Pareto point is of potential interest, the target is to capture the 
whole Pareto front. 

There are several methods for solving multi-objective optimization prob- 
lems; for instance, they can be solved through Evolutionary Algorithms 
(EAs), whose main advantage is that they can find multiple Pareto optimal 
solutions in a single run \W& H] • Among evolutionary methods to tackle 
multi-objective optimization a classical approach entails a form of scalar- 
ization of the criteria vector. Repeated applications of these methods are 
performed to achieve an estimation of the Pareto front. The aggregate objec- 
tive function methods transform a multi-criteria optimization problem into 
a scalar problem using free parameters to be set; for every set of parameter 
values, the scalar optimization problem is solved to seek a Pareto solution. 
Hence, the original problem JTJ is transformed as follows: 

mmG(h(x),...J m (x)), (2) 

xex 

with G : 2 C K m -► R. 

Konak et al. in [T5j provide a tutorial on adopting genetic algorithms to 
solve multi-objective problems, discussing also aggregation approaches based 
on weighted sum of the objective functions. Despite such a method is a quite 
common formulation for a multi-objective optimization problem, there are 
several issues which deserve further investigation; a recent work of Marler 
and Arora [33] proposes further insight on the weighted aggregation method, 
focusing on the choice and significance of the weights according to different 
criteria. Besides, [5] investigates the sensitivity of this class of methods to the 
weights chosen in the scalarization. 

The fundamental issue of a scalarization approach is determining whether 
the transformed problem and the original one are equivalent. In order to 
provide the decision maker with the chance to choose among all optimal 
points, an aggregate function should be able to capture any existing Pareto 
solution. It is possible to prove that any Pareto optimal point can be captured 
if there is an appropriate aggregate function [24j . where a point x is called 
capturable if it is a local optimum of the scalarized problem. Therefore, a main 
issue of this approach is the determination of an appropriate structure for the 
aggregate function, able to cover the Pareto front (as much as possible) by 
properly setting/varying the parameter values. A possible way to deal with 
this aspect is based on the dynamic variation of the aggregation function, 
through properly modifying its parameters during the optimization process. 
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This allows to capture multiple Pareto solutions, reducing at the same time 
the sensitivity issues related to a fixed choice of the weights [5] . 

In this work, an experimental comparative study of different Dynamic 
Objectives Aggregation Methods (DOAMs) in the context of evolutionary 
optimization algorithms is proposed; in particular, we investigate the perfor- 
mance of chaotic variations of the parameters (weights), based on promising 
results from recent research works [6 , 2SJ O [5] on the integration of chaotic 
maps in optimization algorithms. The study is conducted on a set of bench- 
mark problems from the literature. Section[2]presents methods based on both 
weighted sum aggregations, and curvature variations. In Sect. [3] the experi- 
mental setting is described. In Sect. 14, the analysis of results is reported, and 
some conclusions are drawn. 



2 Dynamic Objectives Aggregation Methods 

The aim of this section is the introduction of the evolutionary dynamic ob- 
jectives aggregation methods to solve multi-objective optimization problems. 
Since an aggregate function maps the feasible criteria region into a one- 
dimensional value space, the aggregation method transforms a vector op- 
timization problem into a scalar problem. 

The most common and widely used aggregate function is the weighted sum 
of objectives. Although it has been shown that the weighted sum aggregate 
function is unable to deal with multi-objective optimization problems with 
a concave Pareto front, in (T2J [TS] it is investigated the possibility to cap- 
ture also the points on concave Pareto front by using a dynamic weighted 
aggregation combined with evolution strategies. 

In [16] the phenomenon of global convexity is introduced in order to explain 
the potential success of dynamic weighted aggregation. However, no analyti- 
cal characterization is given to identify a global convex problem, therefore the 
discussion is based on an observed behaviour rather than theoretical analysis. 

From an implementation point of view, classical methods that scalarize 
multiple objectives perform repeated runs in order to achieve a set of non- 
dominated solutions. Dynamic weighted aggregation, instead, provides an 
entire front of non-dominated solutions in a single run. At this aim, these 
procedures generally use an archive to store the non-dominated solutions 
obtained during the search process. 

Empirical results in the literature show that the evolutionary dynamic 
weighted sum is able to provide a set of non-dominated solutions in one 
run of the evolutionary algorithm also capturing in some cases the points 
on concave parts of the Pareto frontier [T3J [151 US]. Besides that, a method 
based on the increase of the aggregate function curvature, obtained varying 
the exponents of the objective function in the aggregation, seems to be able 
to capture the points on concave regions of the front where the plain weighted 
sum fails. 



88 G. Dellino, M. Fedele, and C. Meloni 

The rationale behind an integration of the two methods can be summarized 
observing that by dinamically varying the curvature it may be possible to 
reach the concave part of the front and by dynamically varying the weights 
it may be possible to move close to the concave Pareto frontier. 

The remainder of the section is devoted to introduce the algorithmic ap- 
proaches investigated in our computational study: i) dynamic weighted sum 
methods, and ii) dynamic curvature variations methods. 

2.1 Dynamic Weighted Sum Methods 

The most widely used aggregate function is the weighted sum; the corre- 
sponding aggregate optimization problem can be stated as: 

m 

miny^Wifi(x) , (3) 

X * — » 

i=\ 

where u>i (i = 1, . . . , m) are non- negative weights for the corresponding ob- 
jective functions fi and X}j=i w i = •*- 

For every choice of the weights vector w, the problem ((3| yields an optimal 
Pareto point. It is well-known that a weakness of this aggregate function is the 
failure to capture the points on a concave Pareto fronts. In fact, it is possible 
to prove that every point captured by Yn=i w ifi IS m a convex region of the 
non-dominated frontier. 

In |14j the dynamic weighted aggregation method combined with evolution 
strategies has been studied and it has been shown that this method is able 
to capture the entire Pareto frontier reaching the points in concave regions 
as well. This procedure is based on the dynamic aggregation approach. 

While conventionally the scalarization function weights are fixed during 
optimization, the main idea on which the method is based is that the weights 
systematically change during evolution; so the function to be minimized dy- 
namically changes. In this way the optimizer moves close to the frontier, 
once it achieves a non-dominated solution. During the evolution of the algo- 
rithm, the population can intersect the Pareto set, and therefore an archive is 
needed to record all the Pareto solutions encountered. Although it has been 
extensively shown that the conventional weighted sum is unable to provide 
the Pareto solutions on concave regions, the dynamic weighted sum method 
succeeds in obtaining non-dominated solutions in concave regions as well, 
traversing the frontier dynamically. Empirical results highlight the impor- 
tant role of law of the varying weights P3J [T3] . 

Several ways of changing weights have been proposed: randomly, switching 
between and 1, and periodically. In the first case, the weights are gener- 
ated from a uniform random distribution changing in each generation. The 
second way of varying the weights is realized by switching them from zero 
to one abruptly and viceversa. Literature results suggest that the weights 
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should vary gradually and periodically. In particular, a gradual and continu- 
ous change is needed to keep the points on a convex front: an abrupt switch of 
the search direction does not allow the optimizer to move close to the front, 
storing non-dominated points. 

Since the incorporation of chaos in population-based optimization algo- 
rithms has been shown to possibly enhance their searching ability [S] [T31 [H] , 
this study proposes to introduce and evaluate also the chaotic rules in the 
dynamic weights generation as well. 

2.2 Dynamic Curvature Variation Methods 

In order to overcome the drawbacks of the weighted sum scalarization func- 
tion, several aggregate functions have been introduced in the literature. In 
particular, to enhance the capability of objective functions to capture also 
the points on a concave Pareto front, in [23] it is suggested to increase the 
curvature of the aggregate function. A way of varying the curvature can be 
easily obtained applying the so called i-th power transformation to the objec- 
tive function. The corresponding scalar optimization problem can be stated 
as follows: 



mm 



i=i 



WiifiW , (4) 



where t is a positive real number. It is found that varying all the free param- 
eters (i.e. weights and exponents), it is possible to reach all the points on the 
Pareto frontier. This aggregate function is also investigated in [T51[2D], where 
it has been proved that applying the i-th power to the objective functions the 
convexification of non-dominated frontier can be achieved in an appropriate 
equivalent objectives space. The main problem is again the choice of a func- 
tion structure enabling to provide all the Pareto solutions for some values 
of the parameters used in aggregate function. Assuming that the aggregate 
objective function and the Pareto frontier satisfy certain differentiability re- 
quirements, the necessary and sufficient condition for an admissible aggregate 
objective function to capture the Pareto points has been developed by Messac 
et al. [53] . Although these conditions are inapplicable if the Pareto frontier is 
not known — as it is in real applications — Messac et al. suggested the use 
of an aggregate function (J3J) whose curvature can be increased by setting free 
parameters with the aim to enhance the capability of objective functions to 
capture also the points on concave Pareto front. This i-th power approach is 
also investigated in [JJ5] . For sufficiently large values of t, the efficient frontier 
in the [/*,..., /£J space is guaranteed to be convex under certain conditions. 
Therefore, the weighted sum of the i-th power of the objectives is able to solve 
the problem in the [/*, ••.,/„] space. In [TT] the properties of the weighted 
i-th power approach are summarized: 
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i) the optimal solutions of the i-th power problem ^ are efficient solutions 

of the multi-objective problem ((T|); 
ii) for every efficient solution of the problem (TTJ there exists at > such 

that for all t > i the t-th power aggregate function in (|4j) captures that 

solution. 

This result guarantees the existence of a i-th power aggregate function that 
is able to capture the whole Pareto front. Therefore this is an important the- 
oretical support for our work in which different rules to dynamically change 
the values of t are considered in addition to those concerning the weights Wi. 



3 Computational Experiments 

In this section the experimental setting is illustrated. The evolutionary al- 
gorithms involved in the test and their configurations are described in Sub- 
sect. 13.11 Subsection 13.21 reports on the different DOAMs considered in the 
experiments. In order to evaluate and compare the effectiveness of the pro- 
posed methods, a suite of test problems is employed as will be described in 
Subsect. 



3.1 Evolutionary Algorithms and Their 
Configurations 

To test the methods proposed in this work, the standard genetic algorithm 
included in the Matlab's Genetic Algorithm and Direct Search Toolbox [25] 
has been used. This algorithm enables to solve single-objective optimization 
problems and can be easily adapted to work with dynamic objectives aggre- 
gation. Some parameters values need to be specified, before the algorithm 
execution [29 . In our experiments, we used default (and commonly used) 
values; in particular, we adopted a stochastic uniform selection operator, a 
scattered crossover function with probability 0.7 and a Gaussian mutation 
function with probability 0.3; the number of best individuals that will sur- 
vive to the next population has been fixed to 2 and the stopping criterion 
is based on the maximum number of generations to be produced. A more 
refined tuning will be needed when applying these methods to specific prac- 
tical applications or for a detailed study on the sensitivity of the optimization 
algorithm to these parameters [T] . 

In this work, several alternative settings have been considered for the ge- 
netic algorithm, by varying the population size (in the set {25, 50, 100}) and 
the number of performed iterations (100, 500 or 1000); thus, 9 overall different 
configurations of the genetic algorithm are used. 

The considered DOAMs require to dynamically solve single-objective op- 
timization problems, given by the dynamic weighted sum of the objectives 
we are really interested in. In order to keep all the non-dominated solutions 
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obtained during the optimization process, an archive is needed. At each step 
a dominance analysis on the offspring population is conducted, with respect 
to the functions composing the aggregate objective, thus obtaining a set of 
individuals that are candidated to enter the archive. 

The archive is composed of non-dominated solutions collected as the re- 
sult of the optimization search, that has taken account not only of the ag- 
gregate function but also of its single components. It is important to update 
this archive, removing all dominated solutions. The update frequency of the 
archive (denoted by f up d) is set by the user before the algorithm execution. 

The archive structure A has a maximum size S, but at any iteration it 
may contain a number s = \A\ < S of elements. The use of the archive 
requires a domination analysis, a capacity control and a crowding analysis 
for the elements that are proposed to be enclosed in it. At each iteration the 
evolutionary optimizer proposes to the archive the non-dominated elements 
contained in its current population; they enter the archive according to a 
specific acceptance criterion. More specifically, the following operations are 
performed: a domination analysis ensures that the candidate elements are 
entitled to enter the archive. For each candidate two situations are possible: 
(a) it dominates some of the elements in A, so it replaces the dominated 
solutions; (b) it is non-dominated with respect to all the solutions in A, so it 
could be added to the archive. However, this requires a preliminary capacity 
control, in order to guarantee that the maximum size S is not exceeded. In 
case the archive size reaches the maximum value S, a crowding analysis is 
performed, estimating the density of solutions in the neighborhood of each 
clement of the archive. In particular, a crowding distance value is assigned 
to the individuals belonging to the archive, calculated as follows: for each 
objective function, the solutions in the archive are sorted in ascending order 
of magnitude. Then an infinite distance value is assigned to solutions having 
the smallest and highest objective function values; for solutions other than 
the boundary ones, the distance value is given by 

d k = d k + fi{Xk+l) - fi{Xk - l) , (5) 

J i J i 

where d k is the crowding distance associated to the fc-th clement (fc = 
1, . . . ,\A\), whose initial value is set to 0, fi(-) is the value of the i-th objective 
function of the element specified in brackets, x k +i and Xk-i are the succes- 
sor and the predecessor of the fc-th individual. According to the operated 
sorting, f^ nax and f™ m are the maximum and minimum value of the i-th 
objective function. The calculation of d k is performed with all the objective 
functions, finally obtaining the overall crowding distance for each individual. 
Thus solutions located in lesser crowded regions are preferred. 

Figure 1 outlines how the algorithm works; Figure 2 illustrates the accep- 
tance rule of new elements in the archive. 



92 G. Dellino, M. Fedele, and C. Meloni 



Optimization Process 

Input: Parameters set by the user; 

Initialization: 

Generate a random initial population 
begin 

while the stopping criteria is not satisfied 
Evaluate all individuals in the population 
Rank the population, according to their fitness values 
for each non-dominated individual e in the population 

Archive acceptance rule(e) 
if generation mod f up d = then 

Update the archive, removing dominated elements 
Select parents, based on their fitness 

Create a new population, applying elitism, crossover and mutation operators 
Update the archive, removing dominated elements 
end 



Fig. 1 Genetic Algorithm Scheme. 



Archive acceptance rule 

Input: archive A, archive size S, candidate element e; 

begin 

if e dominates any element in A then 

delete all dominated elements of A and include e in A 
else if no element of A dominates e and \A\ < S then 

include e in A 
else if no element of A dominates e and \A\ = S then 

if the crowding distance of e is better than that of x 6 A then 
delete x and include e in A 
end 



Fig. 2 The acceptance criterion of new elements in the archive. 

A well-known multi-objective genetic algorithm (MOGA), NSGA-II - ex- 
tensively described by Deb et al. in [10] - has also been used, aiming to 
compare the solutions obtained with the proposed method with those pro- 
vided by a native MOGA. The algorithm configurations described before 
have also been applied to NSGA-II, except for the dominance management 
which is implicitly guaranteed by the algorithm itself; NSGA-II also includes 
a crowding analysis similar to the one we adopt (the details of NSGA-II can 
be found in [ID]). 
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3.2 The Set of DOAMs 

Several DOAMs have been used in the campaign of experiments conducted 
in this work, involving both the variation of the weights only, as in ((3J), and 
the combined variation of the weights and the exponents in ((H) . For the sake 
of simplicity, in order to illustrate the methods, a bi-objective problem is 
considered. In this case, the aggregate function corresponding to the fc-th 
generation can be stated as follows: 

G(x, k) = w x {k)f\{x) + w 2 (k)f 2 (x) , (6) 

where the expressions of w\ and w 2 and the value of t depend on the adopted 
variation law; clearly, for t — 1 the simpler weighted sum aggregated function 
is obtained. The weights Wi can be dynamically modified according to a rule 
R(k) described by a specific function of k: 

wi(k) = R(k), w 2 (k) = 1 - wi(k) . (7) 

The weighted aggregation methods can easily be extended to three-objective 
problems |17| ; the weights can be generated in the following way: 

wi(k) = R(k), w 2 {k) = (1- wi(k))R(k), w 3 (k) = 1 - wi(k) - w 2 {k) . 

We consider different rule functions, namely: one, switch, sin, triangle, rand, 
chaos. 

The first refers to the case of fixed unitary weights (i.e. the aggregate 
function is simply given by the sum of the objectives). In the second case w\ 
periodically changes from to 1, with a given period T = 100 (in terms of 
the number of algorithm's generations). Similarly, a periodical changing of 
the weights can be obtained also according to a sin or triangle wave in the 
successive adopted rules; the sinusoidal rule is the following: 

R(k) = \sm(2wk/F)\, (8) 

where F is the frequency in terms of algorithm's generations. In our com- 
putational experiments, this parameter has been set to F — 200 following 
|14| showing that such a value leads to algorithms with better convergence 
properties. The rand rule, at each iteration k, generates a random value in 
(0,1) for wi. The last rule applies a chaotic variation law to the weights. 
A logistic equation - which is extensively used to describe a chaotic system 
[s2 O HI] - is employed as follows: 

wi(* + l) = Aiu>i(*)(l-«;i(fc)), w 2 {k) = 1 - wi(k) ; (9) 

where fi = 4 and itfi(O) = 0.2027. The previous deterministic equation (choos- 
ing tui(0) 6 (0, 1)\{0, 0.25, 0.5, 0.75, 1}) shows chaotic behaviour, exhibiting a 
sensitive dependence on initial conditions, which is the basic characteristic of 
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chaos. A little difference in the initial value of the chaotic variable would result 
in a considerable difference in its long-time behaviour. In general, this chaotic 
variable has special characteristics, such as ergodicity, pseudo-randomness 
and irregularity [H 03 [21] . Clearly, some other well-known chaotic maps could 
also be employed instead of the logistic one to generate the weights in the 
aggregate objective functions [5]. 

In order to let the curvature of the aggregate function vary during the 
evolution process four possible strategies are proposed for the variation of 
the exponent. In all the cases considered in the following, the exponent value 
ranges between t = 1 and t = 4, retaining that greater values of t would not 
provide further improvements in the optimization results achieved so far. 

A first scheme (one) considers only fixed unitary exponents. The second 
scheme (step) establishes to increment the exponent value every JV/4 itera- 
tions, N being the maximum number of generations that can be produced. 
An adaptive scheme (adapt) has also been tested, according to which the 
exponent value is incremented when there is no improvement in the opti- 
mization process for a given number of iterations, which has been fixed to 
A = 0.05N. 

According to both these strategies, the exponent value is always a posi- 
tive integer number. The last strategy (cont) considered in this work let the 
exponent range among the (positive) real numbers; i.e., the interval (0, N) 
has been mapped into the interval (1,4) such that the exponent t can vary 
continuously in this range. 

Combining the aforementioned weights-exponents strategies, 24 different 
algorithms are obtained. Hereinafter, each of them will be denoted indicating 
the strategies as an ordered pair (e.g. chaos-step represents the strategy with 
the chaotic rule for the weights and the step rule for exponents, respectively) 
and they will be labeled as reported in Table [TJ 

Table 1 The considered DOAMs and the corresponding weights-exponents 
strategies 



Algorithm Weights-Exponents Rule 


Algorithm 


Weights-Exponents Rule 


DOAMl 


chaos-one 


DOAM13 


chaos- adapt 


DOAM2 


one-one 


DOAM14 


one- adapt 


DOAM3 


rand-one 


DOAM15 


rand-adapt 


DOAM4 


switch-one 


DOAM16 


switch- adapt 


DOAM5 


sin- one 


DOAM17 


sin-adapt 


DOAM6 


triangle-one 


DOAM18 


triangle- adapt 


DOAM7 


chaos-step 


DOAM19 


chaos-cont 


DOAM8 


one-step 


DOAM20 


one-cont 


DOAM9 


rand-step 


DOAM21 


rand-cont 


DOAM10 


switch-step 


DOAM22 


switch-cont 


DOAM11 


sin-step 


DOAM23 


sin- cont 


DOAM12 


triangle-step 


DOAM24 


triangle-cont 



Dynamic Objectives Aggregation Methods 95 

3.3 Test Problems 

The computational test of the methods has been conducted on a set of bench- 
mark problems, characterized by different specific features in the Pareto front, 
so that the general results obtained would not depend on the particular test 
problem chosen. Problems P1-P7 are discussed by Jin et al. in [HI [T51 fT5] : 
in P2-P5 it is assumed that Xi £ [0, 1] for all i = 1, . . . ,n; while in Pi, Pg, 
and P7 there are no restrictions on the range of the decision variables. The 
problem Pi has the following objective functions: 

_. n 1 n 

h = -J2 x i> / 2 = -]>>»- 2 ) 2 ; (10) 

n * — ' n * — ' 

and produces a uniform Pareto front. P2 is described by: 

h = x i 

9 <A 

g(x 2 ,...,x n ) = 1 + 7>£i, (11) 

n — 1 * — ' 

i=2 



h = .9 x (1 - y/h/g) . 

having a convex Pareto front. Because of the interest in studying problems 
showing non-convex or discontinuous Pareto front, some instances belonging 
to this class have been considered. The following problem, P3, has a concave 
Pareto front and is defined as follows: 

/1 = %i , 

9 <A 

g(x2,...,x n ) = H TZ^ X *> ( 12 ) 

n — 1 * — ' 

i=2 

h = 9 x (1 - {h/gf) . 

The fourth problem, P4, has been obtained through combining — in some 
sense — P2 and P3: 

/1 = xi , 

9 " 

g(x 2 ,...,x n ) = H 7^^ 1 ( 13 ) 

ri — 1 * — ' 



/ 2 =.gx(l-^y7 1 7g-(/ 1 /. 9 ) 4 ). 

Its Pareto front is neither purely convex nor purely concave. The following 
problem, P5, is characterized by a Pareto front consisting of a number of 
separated convex parts. 
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fl = x i > 

9 ™ 

g(x 2 ,...,x n ) = l+ -Vii, (14) 

n — 1 ^^ 

/ 2 = .9 x (1 - yfhfg- (h/g) onClOTr/x)) . 

Problem Pg is defined through the following equations: 

2' 



(15) 



showing a concave Pareto front. Another problem, P7, is taken from [IB] 
extending one of the test function proposed in the literature [55j to the two- 
dimensional case: 

/1 = exp(— xi) + 1.4 exp(— x\) + exp(— x 2 ) + 1.4 exp(— x\) , 

f 2 = exp(zi) + 1.4 exp(-xi) + exp(a;2) + 1.4 exp(-iE2) ■ (16) 

The resulting Pareto front is continuous and non-convex. Even if the problem 
is considered an easy task for evolutionary optimizers |16) . the region that 
defines the Pareto front in the parameter space is disconnected; so, it could be 
an interesting problem to be studied. In the following, two other benchmark 
problems, P$-P§, are considered, because of the particular shape of their 
feasible region and/or Pareto front; these problems are described in [27] . 
Problem Pg is referred to as the TNK problem. The objectives are very 
simple, and defined by 

/1 =xi, h = x% , (17) 

where 

xi e [0,7r], x 2 e [0, 7r] . 

The constraints are 

g\ = —xi — x\ + 1 + 0.1 cos(16arctan(;E2/£i)) < , 

.92 = (xi - 0.5) 2 + (a? 2 - 0.5) 2 < 0.5 . (18) 

Problem Pg is the so-called Poloni's test problem: the objective functions are 
defined as 

h = l + {a-bf + {c-df , h = {xx + 3) 2 + (32 + l) 2 , (19) 
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where the parameters introduced in the expression of f± are: 

a = 0.5sin(l) - 2.0cos(l) + 1.0sin(2) - 1.5cos(2) , 

b = 0.5sin(xi) — 2.0cos(xi) + 1.0 sin^) — 1.5cos(x2) , 

c = 1.5 sin(l) - 1.0 cos(l) + 2.0 sin(2) - 0.5 cos(2) , 

d = 1.5sin(xi) — l.Ocos(xi) + 2.0 sin^) — 0.5cos(x2) . 

The variables ranges are: x\ G [— 7r,7r] and X2 £ [— 7r, tt]. 

4 Discussion of Results 

Appropriate metrics must be selected in order to evaluate the behaviour of 
the considered algorithms. The literature offers different metrics to measure 
the performance of algorithms for multi-criteria optimization problems. Nev- 
ertheless, no single metric is able to assess the total algorithmic performance. 
The metrics adopted in this study are reported below. Clearly, these met- 
rics should not be considered as a complete list of all possible metrics. For 
instance, in our computational experiments, we do not consider particular 
temporal metrics, limiting our analysis only to the computation times re- 
quired by the algorithm. Although used with cases with few objectives, the 
considered metrics can also be applied when a larger number of objectives is 
considered. 

We are interested in measuring how far the non-dominated solutions ob- 
tained by the algorithm, i.e. the solution front (SF), are from the ideal point 
(IP) . An ideal point (IP) is defined as a point characterized by the best val- 
ues for each objective. We also use the concept of the nadir point (NP) which 
is defined as a point characterized by the worst values for each objective. To 
this aim, the adopted measure is given by an hyperarea or hypervolume ratio. 
This metric requires the knowledge (i.e. the computation) of an ideal point 
(IP) and a nadir point (NP) for the problem. The ideal point and the nadir 
point can be viewed as vertices of a regular polytope defining a hypervolume 
(A t ), i.e. the total area, in the space of the objectives (a rectangle in the two- 
objective case). The hyperarea ratio (HR) performance measure is defined as 

HR = Ai/ At , (20) 

where Ad indicates the dominated area between the nadir point NP and 
the solution front SF, as proposed by Fleischer [12,. A large value of HR 
is expected from a good algorithm. The hyperarea ratio is a good perfor- 
mance indicator, but it does not take into account how the efficient points 
are distributed on the estimated solution front. 

Therefore as secondary indicators, we report the number of non- dominated 
elements (ND) and the spacing (S). The last one is a metric measuring the 
spread (distribution) of vectors throughout the solution front (SF) . We refer 
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to the definition reported in [3T] for two objective functions: the spacing 
measures the range variance of neighboring vectors in SF: 



S 



\ 



N 

N -~ 



1 N 

^(d-ck) 2 , (21) 



where d l = minj(\f{(x) - f((x)\ + |/vj(«) - f 2 {x)\), i,j = l,...,N;dis the 
mean of all di, and N is the number of vectors in SF. A value of zero for 
this metric indicates that all non-dominated solutions are equally spaced. It 
is worth observing that no single metric can completely capture total multi- 
objective evolutionary algorithm performance |31| . 

For each experiment, three different runs have been executed initializing 
algorithms with random populations. Tables [2] and [3] contain the average 
results on 3 runs for each of the 9 algorithms configurations described in 
Subsect. 13.11 so each entry is the average over 27 experimental values. 

Table 2 Average values of HR (in %), ND, and S achieved by the algorithms for 
problems P1-P5 



Algorithms 


HR 


Pi 

ND S 


HR 


P2 

ND 


S 


P3 
HR ND S 


Pi 

HR ND S 


P5 
HR ND S 


DOAMl 


74 


25.37 0.36 


90 


10.67 0.42 


85 


6.26 0.17 


81 


7.00 0.29 


86 


6.04 0.21 


DO AM 2 


74 


24.85 0.38 


91 


7.59 


0.45 


83 


4.33 0.29 


81 


7.04 0.25 


83 


5.19 0.40 


DOAM3 


75 


25.41 0.29 


91 


8.52 


0.32 


84 


4.96 0.35 


81 


6.33 0.29 


84 


5.26 0.56 


DOAM4 


73 


23.52 0.53 


88 


3.48 


0.18 


82 


2.22 0.13 


81 


5.78 0.29 


83 


3.78 0.42 


DO AM 5 


75 


24.93 0.31 


91 


8.07 


0.24 


85 


4.44 0.18 


81 


6.37 0.25 


86 


5.07 0.33 


DOAM6 


76 


24.15 0.33 


90 


6.11 


0.37 


85 


5.22 0.18 


81 


6.74 0.27 


85 


5.74 0.36 


DO AM 7 


75 


25.07 0.36 


90 


8.04 


0.30 


84 


5.11 0.32 


81 


6.41 0.25 


85 


6.44 0.36 


DOAM8 


75 


25.00 0.38 


90 


7.78 


0.28 


84 


4.04 0.30 


81 


7.48 0.31 


83 


5.41 0.55 


DOAM9 


76 


24.70 0.31 


91 


6.93 


0.40 


84 


4.33 0.29 


81 


6.67 0.27 


85 


5.22 0.36 


DOAM10 


75 


23.52 0.46 


89 


3.67 


0.16 


82 


2.59 0.23 


80 


5.48 0.35 


84 


3.85 0.34 


DOAM11 


75 


24.63 0.32 


91 


7.22 


0.28 


84 


4.11 0.17 


81 


7.15 0.29 


85 


5.44 0.27 


DOAM12 


75 


25.44 0.28 


91 


8.33 


0.31 


84 


4.67 0.31 


81 


7.52 0.27 


85 


5.33 0.32 


DOAM13 


74 


24.85 0.32 


91 


8.96 


0.20 


84 


4.89 0.21 


81 


7.15 0.26 


85 


5.78 0.28 


DOAM14 


76 


25.93 0.32 


91 


7.41 


0.25 


84 


3.89 0.25 


81 


7.44 0.28 


83 


5.33 0.60 


DOAM15 


73 


25.33 0.31 


90 


6.07 


0.30 


85 


4.70 0.14 


81 


6.78 0.31 


83 


6.15 0.42 


DOAM16 


73 


24.26 0.59 


88 


3.04 


0.13 


83 


3.00 0.18 


80 


4.74 0.37 


84 


4.07 0.32 


DOAM17 


75 


25.30 0.27 


91 


6.67 


0.13 


85 


4.85 0.14 


81 


6.74 0.32 


86 


5.48 0.41 


DOAM18 


76 


25.15 0.38 


91 


8.52 


0.25 


84 


4.56 0.32 


82 


7.44 0.23 


85 


6.26 0.31 


DOAM19 


75 


24.93 0.30 


90 


8.52 


0.25 


84 


5.22 0.41 


81 


7.11 0.29 


85 


6.63 0.27 


DOAM20 


75 


25.22 0.33 


91 


7.93 


0.23 


84 


5.07 0.37 


81 


7.11 0.31 


85 


5.74 0.40 


DOAM21 


75 


23.89 0.38 


91 


7.07 


0.31 


83 


4.19 0.31 


81 


6.78 0.30 


85 


6.11 0.43 


DOAM22 


74 


23.85 0.48 


89 


3.74 


0.26 


82 


3.04 0.10 


80 


5.48 0.32 


83 


3.96 0.12 


DOAM23 


75 


24.81 0.33 


91 


7.63 


0.36 


84 


5.48 0.58 


81 


7.15 0.25 


85 


5.89 0.32 


DOAM24 


76 


25.59 0.33 


90 


7.30 


0.42 


84 


5.85 0.30 


81 


7.19 0.22 


85 


5.89 0.49 


NSGA-II 


65 


26.52 0.07 


90 


4.19 


0.17 


84 


2.59 0.15 


81 


4.96 0.21 


86 


4.11 0.25 
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The collected results show that methods based on periodical changes of 
weights often achieve relatively good performances with respect to other 
DOAMs as well as to a state-of-the-art (native) multi-objective optimizer. 
This behaviour seems to be reinforced by the use of exponents variation. 
Nonetheless, it is noticeable the competitive performance of the DOAMs 
based on a chaos rule in terms of HR, ND, and S. On the other hand, 
strategies based on the switch rule (no matter what strategy is adopted for 
exponents) almost always give relatively modest results. 

Although each average is composed of a large number of data points, it 
is necessary to carry out a statistical analysis to assess if the observed dif- 
ferences in the average values are indeed statistically significant. Even if in 
the optimization literature, parametric test methods are not often adopted 
due to the assumptions that the data have to satisfy and nonparametric test 
methods are in general preferred because they are distribution-free, the for- 
mer ones result more powerful allowing for a much deeper analysis of data. In 
fact, in non-parametric testing a lot of information is lost because the data 

Table 3 Average values of HR (in %), ND, and S achieved by the algorithms for 
problems P6-P9 

~ ~ P^ P~r Pi P9 

Algorithms HR Np S H R ND S HR ND S HR ND S 

16.79 0.14 100 328.19 3.41 23 4.70 0.18 91 224.79 0.21 

15.29 0.14 100 152.93 1.5c56 21 4.89 0.12 89 126.25 0.63 

16.13 0.18 100 310.37 3.68 22 4.63 0.14 91 212.92 0.27 

14.96 0.14 100 474.59 1.7ell 21 4.48 0.18 82 36.63 1.56 

19.17 0.11100 371.56 8.98 22 4.67 0.16 91 160.58 1.18 

20.46 0.10 100 327.00 4.82 21 4.48 0.16 91 166.29 1.29 

17.46 0.16 100 329.52 5.31 21 4.11 0.14 91 244.50 0.23 

15.08 0.18 100 141.93 3.54 22 5.00 0.13 91 158.46 0.27 

17.04 0.12 100 315.37 3.78 21 4.37 0.16 91 231.00 0.26 

15.92 0.12 100 475.74 1.8el0 22 4.59 0.15 83 38.71 1.60 

20.83 0.15 100 373.44 9.37 21 4.30 0.18 91 176.75 0.68 

20.21 0.11 100 325.22 3.70 22 4.15 0.12 91 184.75 0.60 

17.96 0.11 100 339.44 3.90 22 4.48 0.14 91 243.42 0.41 

16.75 0.14 100 144.70 4.08 22 4.78 0.13 91 142.42 0.27 

16.79 0.12 100 319.78 3.92 21 4.15 0.12 91 218.13 0.24 

15.63 0.14 100 478.63 1.3el2 21 4.11 0.16 85 39.46 0.96 

21.21 0.08 100 375.22 11.70 20 4.30 0.14 91 171.42 0.89 

19.13 0.14 100 332.15 4.19 22 4.59 0.19 91 196.75 0.91 

17.08 0.13 100 287.00 4.55 21 4.15 0.12 91 210.83 0.26 

15.54 0.20 100 93.52 3.27 22 4.48 0.15 91 148.50 0.29 

17.17 0.14 100 256.85 4.45 23 4.56 0.15 91 184.96 0.27 

16.38 0.13 100 482.89 2.9el0 21 4.30 0.11 80 32.25 1.17 

19.25 0.10 100 290.78 5.85 22 4.48 0.17 91 162.50 1.08 

20.17 0.11 100 258.74 3.49 22 4.33 0.13 91 171.58 1.20 

NSGA-11 18 16.92 0.03 100 411.33 1.5ell 16 2.48 0.12 91 440.42 0.16 
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DOAM1 - 
DOAM2 - 
DOAM3 - 
DOAM4 - 
DOAM5 - 
DOAM6 - 
DOAM7 - 
DOAM8 - 
DOAM9 - 
DOAM10- 
DOAM11 - 
DOAM12- 
DOAM13- 
DOAM14- 
DOAM15- 
DOAM16- 
DOAM17- 
DOAM18- 
DOAM19- 
DOAM20 - 
DOAM21 - 
DOAM22 - 
DOAM23 - 
DOAM24 - 
NSGA-II - 



Fig. 3 Means plot and Tukey's HSD confidence intervals (a — 0.05) resulting from 
the Rank-based Friedman analysis on HR 

have to be ranked and the differences in the values are transformed into a 
rank value; therefore, besides the non-parametric rank-based Friedman's test, 
we consider parametric ANOVA analysis [7] |Hl US] ■ 

In these analyses we assume only HR as response variable. 

Figures[3jand|4]show the means plot in the Friedman and ANOVA analysis, 
respectively. The analyses are conducted on the algorithm factor considering 
24 different DOAMs and NSGA-II with three replicates for each experiment 
which is characterized by the problem and by the algorithm configuration 
(a total of 81 combinations are considered). Thus, in our Friedman analysis, 
for every experiment 75 ranks are obtained, assigning a larger rank to better 
results. For both tests we use a 95% confidence interval and adopt the Tukey's 
HSD intervals [30] US] • As it can be seen in Figures [3J and QJ the non- 
parametric test is less powerful neglecting the differences in the response 
variables and presenting much wider confidence intervals. In the reported 
plots, overlapping confidence intervals indicate a non-statistically significant 
difference on the average performance of the respective algorithms. 

These analyses clearly confirm the negative assessment on switch strategies 
and the promising behaviour of chaos-based DOAMs which are often the best 
strategy (e.g. see DOAM13 in Figure [3]), providing good performance based 
on the three metrics considered in Tables [2] and |3j and being non-dominated 
with respect to the other strategies (see Figures [3] and |4| . This result seems 
to support the increasing research interest in the introduction of some form 
of chaotic behaviour in stochastic optimizers [3J [5] [T3J [55] . 

These experimental results are of interest also in different contexts such 
as: the development of multi-objective optimization algorithms starting from 
well-established evolutionary single-objective optimizers; the design of com- 
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Fig. 4 Means plot and Tukey's HSD confidence intervals (a = 0.05)resulting from 
the ANOVA analysis on HR 

pact (and fast) local search procedures; surrogate-based optimization; land- 
scape approximation of costly functions. Moreover, the case of bi-objective 
optimization could be of interest in Mean- Variance optimization problems 
often used to model either the risk preferences of the decision maker or ro- 
bustness requirements. 

The observations based on the encouraging results from the conducted 
experiments indicate that different aspects deserve further research efforts 
including the extension to the experimental campaign on a wider set of prob- 
lems (also from real applications); the possible use of different quality indi- 
cators; the comparison of DOAMs with other state-of-the-art MOGAs; the 
investigation of the effects of the interactions of weights and exponents based 
rules; and the consideration of other chaotic DOAMs in order to deeply in- 
vestigate their effectiveness. 
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Summary. Network-on-chip (NoC) are considered the next generation of commu- 
nication infrastructure, which will be omnipresent in most of industry, office and 
personal electronic systems. In the platform-based methodology, an application is 
implemented by a set of collaborating intellectual properties (IPs) blocks. In this 
paper, we use multi-objective evolutionary optimization to address the problem of 
mapping topologically pre-selected sets IPs, which constitute the set of optimal so- 
lutions that were found for the IP assignment problem, on the tiles of a mesh-based 
NoC. The IP mapping optimization is driven by the area occupied, execution time 
and power consumption. 



1 Introduction 

As the integration rate of semiconductors increases, more complex cores for system- 
on-chip (SoC) are launched. A simple SoC is formed by homogeneous or heteroge- 
neous independent components while a complex SoC is formed by interconnected 
heterogeneous components. The interconnection and communication of these com- 
ponents form a network-on-chip (NoC). A NoC is similar to a general network but 
with limited resources, area and power. Each component of a NoC is designed as an 
intellectual property (IP) block. An IP block can be of general or special purpose 
such as processors, memories and DSPs |7|. 

Normally, a NoC is designed to run a specific application. This application, usu- 
ally, consists of a limited number of tasks that are implemented by a set of IP blocks. 
Different applications may have a similar, or even the same, set of tasks. An IP block 
can implement more than a single task of the application. For instance, a processor 
IP block can execute many tasks as a general processor does but a multiplier IP block 
for floating point numbers can only multiply floating point numbers. The number of 
IP blocks designers, as well as the number of available IP blocks, is growing up fast. 
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In order to yield an efficient NoC-based design for a given application, it is nec- 
essary to choose the adequate minimal set of IP blocks. With the increase of IP 
blocks available, this task is becoming harder and harder. Besides IP blocks care- 
fully assignment, it is also necessary to map the blocks onto the NoC available 
infra-structure, which consists of a set of cores communicating through switches. 
A bad mapping can degrade the NoC performance. Different optimization criteria 
can be pursued depending on how much information details is available about the 
application and IP blocks. 

Usually, the application is viewed as a graph of tasks called task graph (TG). The 
IP blocks features can be obtained from their companion documentation. The IP 
assignment and IP mapping are key research problems for efficient NoC-based de- 
signs 1151 . These two problems can be solved using multi-objective optimizations 
in which some of the objectives are conflicting. Because of their nature, both IP 
assignment and mapping are classified as NP-hard problems J6|. Normally, deter- 
ministic techniques are not viable to solve such problems, so we use multi-objective 
evolutionary algorithms (MOEAs) with specific operators and objective functions 
to yield an optimal IP mapping for a prescribed set of IP assignments. These consti- 
tute the set of optimal solutions that were found in the IP assignment satge. For this 
purpose, one needs to select the best minimal set of objectives to be optimized. The 
wrong set of optimized objectives can lead to average instead of best results. Here, 
we assume that the IP assignment has been performed and is available for mapping 
the application TG onto the NoC infrastructure. 

In this paper, we propose a multi-objective evolutionary-based decision support 
system to help NoC designers. For this purpose, we propose a structured repre- 
sentation of the TG and an IP repository that will feed data into the system. We 
use the data available in the Embedded Systems Synthesis benchmarks Suite (E3S) 
l3ll as our IP repository. The E3S is a collection of TGs, representing real appli- 
cations based on embedded processors from the Embedded Microprocessor Bench- 
mark Consortium (EEMBC). It was developed to be used in system-level allocation, 
assignment, and scheduling research. We used two MOEAs: NSGA-II J2) and mi- 
croGA 1 1 1. Both of these algorithms were modified according to some prescribed 
NoC design constraints. 

The rest of the paper is organized as follows: First, in Section |2j we present 
briefly some related research work. Then, in Section [3J we introduce an overview 
of NoC structure. Subsequently, in Section [4] we describe a structured TG and IP 
repository model based on the E3S data. After that, in Section [5j we introduce the 
mapping problem in NoC-based environments. Then, in Section [7] we sketch the 
two MOEAs used in this work, individual representations and objective functions 
for the optimization stage. Later, in Section [9] we show some experimental result. 
Last but not least, in Section \W\ we draw some conclusions and outline new direc- 
tions for future work. 
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2 Related Work 

The problems of mapping IP blocks into a NoC physical structure have been ad- 
dressed in some previous studies with different emphasis. Some of these works did 
not take into account of the multi-objective nature of these problems and adopted a 
single objective optimization approach. 

Hu and Marculescu Q proposed a branch and bound algorithm which automat- 
ically maps IPs/cores into a mesh based NoC architecture that minimizes the total 
amount of consumed power by minimizing the total communication among the used 
cores. Specified constraints through bandwidth reservation were defined to control 
communication limits. 

Lei and Kumar ill 1 proposed a two step genetic algorithm for mapping the TG 
into a mesh based NoC architecture that minimizes the execution time. In the first 
step, they assumed that all communication delays are the same and selected IP 
blocks based on the computation delay imposed by the IPs only. In the second step, 
they used real communication delays to find an optimal binding of each task in the 
TG to specific cores of the NoC. 

Murali and De Micheli |[T31 addressed the problem under the bandwidth con- 
straint with the aim of minimizing communication delay by exploiting the possibil- 
ity of splitting traffic among various paths. Splitting the traffic increases the size of 
the routing component at each node but the authors were not worried about size. 

Zhou et al. ET1 suggested a multi-objective exploration approach, treating the 
mapping problem as a two conflicting objective optimization problem that attempts 
to minimize the average number of hops and achieve a thermal balance. The num- 
ber of hops is incremented every time a data cross a switch before reaching its 
target. They used NSGA ifTHl . multi-objective evolutionary algorithm. They also 
formulated a thermal model to avoid hot spots, which are areas with high computing 
activity. 

Jena and Sharma |9| addressed the problem of topological mapping of IPs/cores 
into a mesh-based NoC in two systematic steps using the NSGA-II (2|. The main ob- 
jective was to obtain a solution that minimizes the energy consumption due to both 
computational and communicational activities and also minimizes the link band- 
width requirement under some prescribed performance constraints. 

As a recent field of research, the available literature about NoC-based design 
optimization is scarce. The aforementioned works represent the state of the art of 
this field. In Q, ITU and llT3l . only one objective was considered and only ifTTI 
used an evolutionary approach. In 1*211 and J9), two objectives were considered and 
both adopted a MOEA to solve the problem. 

3 NoC Internal Structure 

A NoC platform consisting of architecture and design methodology, which scales 
from a few dozens to several hundreds or even thousands of resources ITOl . As 
mentioned before, a resource may be a processor core, DSP core, an FPGA block, 
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a dedicated hardware block, mixed signal block, memory block of any kind such as 
RAM, ROM or CAM or even a combination of these blocks. 

A NoC consists of set of resources (R) and switches (S). Resources and switches 
are connected by channels. The pair (R, S) forms a tile. The simplest way to connect 
the available resources and switches is arranging them as a mesh so these are able to 
communicate with each other by sending messages via an available path. A switch 
is able to buffer and route messages between resources. Each switch is connected to 
up to four other neighboring switches through input and output channels. While a 
channel is sending data another channel can buffer incoming data. Note that energy 
consumption is proportional to the number of message exchanges. Therefore, the 
communication between resources and distance between them must be considered 
during the mapping stage. Fig.Q]shows the architecture of a mesh-based NoC where 
each resource contains one or more IP blocks (RNI for resource network interface, D 
for DSR M for memory, C for cache, P for processor, FP for floating-point unit and 
Re for reconfigurable block). Besides the mesh topology, there are more complex 
topologies like torus, hypercube, 3-stage clos and butterfly lfT4ll , Note that mesh- 
based NoC does not always represent the best topological choice. 




Fig. 1 Mesh-based NoC with 9 resources 



Every resource has an unique identifier and is connected to the network via a 
switch. It communicates with the switch through the available RNI. Thus, any set 
of IP blocks can be plugged into the network if its footprint fits into an available 
resource and if this resource is equipped with an adequate RNI. The NoC imple- 
ments the lower four layers from OSI model |8|: physical, data-link, network and 
transport layers. The RNI implements all four layers and it packs transport layer 
messages into network layer packets; so, the switch-to-switch interface implements 
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only the three lower protocol layers. The basic communication mechanism envi- 
sioned among computing resources is message passing. However, it is possible to 
add additional protocols on top of the transport layer. 

4 Task Graph and IP Repository Models 

In order to formulate the IP mapping problem, it is necessary to introduce a formal 
definition of an application first. An application can be viewed as a set of tasks that 
can be executed sequentially or in parallel. It can be represented by a directed graph 
of tasks, called task graph. 

Definition 1. A Task Graph (TG) G = G(T, D) is a directed graph where each node 
represents a computational module in the application referred to as task ti G T. Each 
directed arc dij G D, between tasks ti and tj, characterizes either data or control 
dependencies. 

Each task ti is annotated with relevant information, such as a unique identifier and 
type of processing element (PE) in the network. Each dij is associated with a value 
V(dij), which represents the volume of bits exchanged during the communication 
between tasks ti and tj . 

Once the IP assignment has been completed, each task is associated with an IP 
identifier. The result of this stage is a graph of IPs representing the PEs responsible 
of executing the application. 

Definition 2. An Application Characterization Graph (APG) G = G(C, A) is a 
directed graph, where each vertex a G C represents a selected IP/core and each 
directed arc a,,j characterizes the communication process from core Ci to core cj. 

Each a,.j can be tagged with IP/application specific information, such as commu- 
nication rate, communication bandwidth or a weight representing communication 
cost. 

A TG is based on application features only while an APG is based on application 
and IP features, providing us with a much more realistic representation of the an 
application in runtime on a NoC. In order to be able to bind application and IP 
features, at least one common feature is required in both of the IP and TG models. 

The E3S (0.9) Benchmark Suite |3j contains the characteristics of 17 embedded 
processors. These processors are characterized by the measured execution times of 
47 different type of tasks, power consumption derived from processor datasheets, 
and additional information, such as die size, price, clock frequency and power con- 
sumption during idle state. In addition, E3S contains task graphs of common tasks 
in auto-industry, networking, telecommunication and office automation. Each one 
of the nodes of these task graphs is associated with a task type. A task type is a 
processor instruction or a set of instructions, e.g., FFT, inverse FFT, floating point 
operation, OSPF/Dijkstra 0, etc. If a given processor is able to execute a given 
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type of instruction, so that processor is a candidate to receive a resource in the NoC 
structure and would be responsible for the execution of one or more tasks. 



4.1 XML Representation 

The E3S Benchmark Suite contains rich data about embedded processors and some 
common applications. TGFF [4], a random TG generator based on E3S processors, 
generates TGs with parallel and sequential tasks, nodes with IP types and other im- 
portant features. Both, E3S and TGFF, are text files. We use XML Schema EU| 
to model the TG and IP repository. At this point, no standard schema for NoC de- 
sign is available, so the XML structure for both representations reflects the features 
available from E3S processors and applications. 

XML is a general-purpose specification for creating custom markup languages 
and we propose to use it as a standard in NoC design research. Its primary purpose 
is to help information systems share structured data and it is designed to be relatively 
human-legible. It is strongly structured and its structure can be controlled by a XSD 
schema. Any XML file based on a schema can be readable for any tool designed for 
that schema. To parse (process) a XML file is much easier than a TXT file. Modern 
programming languages offer APIs that facilitate XML files parsing, while TXT 
files must be read line by line, checking character by character. XML and XML 
schema (XSD) are defined by the World Wide Web Consortium (W3C) (2D|. 

4.2 Task Graph Representation 

Here, we represent TGs using the XML l20ll . A TG is divided in three major ele- 
ments: taskgraph, nodes and edges. The taskgraph element is the TG itself which 
contains nodes and edges. The nodes element includes a node element for each task 
of the TG and the edges element includes an edge element for each arc in the TG. 
Each node has two main attributes: an unique identifier (id) and a task type (type), 
chosen among the 47 different types of tasks present in the E3S. Each edge has four 
main attributes: an unique identifier (id), the id of its source node (src), the id of 
its target node (tgt) and an attribute representing the communication cost imposed 
(cost). Fig.|2]shows a simple TG and its corresponding XML representation. 



4.3 IP Repository Representation 

The IP repository is divided into two major elements: the repository and the ips 
elements. The repository is the IP repository itself. Recall that the repository con- 
tains different non general purpose embedded processors and each processor imple- 
ments up to 47 different types of operations. Not all 47 different types of operations 
are available in all processors. Each type of operation available in each processor 
is represented by an ip element. Each ip is identified by its attribute id, which is 
unique, and by other attributes such as taskType, taskName, taskPower, taskTime, 
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<taskgraph> 










<nodes> 










<node 


id= 


"0" 


type="45" . . . /> 




<node 


id= 


"1" 


type="21" . . . /> 




<node 


id= 


"2" 


type="28" . . . /> 




<node 


id= 


"3" 


type="32" . . . /> 




</nodes> 










<edges> 










<edge 


id= 


"0" 


src="0" tgt="l" cost= 


"5"/> 


<edge 


id= 


"1" 


src="l" tgt="3" cost= 


:"4"/> 


<edge 


id= 


"2" 


src="l" tgt="2" cost= 


:"3"/> 


<edge 


id= 


"3" 


src="2" tgt="3" cost= 


:"2"/> 


</edges> 










</taskgraph> 









Fig. 2 TG XML structure 

processorlD, processorName, processorWidth, processorHeight, processorClock, 
processorldlePower and cost. The common element in TG and IP repository rep- 
resentations is the type attribute. Therefore, this element will be used to bind an 
ip to a node. Fig. [3] shows a simplified XML structure representing the IP repos- 
itory. The repository contains IPs for digital signal processing, matrix operations, 
text processing and image manipulation. 

<?xml version= " 1 . " encoding= "UTF-8 " ?> 



<repository> 












<ips> 












<ip id= 


"10" 


type = 


11 0" 


procID="3 " . . 


./> 


<ip id= 


"23" 


type = 


"38' 


' procID="5" 


. ./> 


<ip id= 


"68" 


type= 


"12' 


' procID="14" 


. . ./> 


<ip id= 


"99" 


type= 


"47' 


' procID="17" 


. . ./> 


</ips> 













</repository> 



Fig. 3 IP repository XML structure 
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These simplified and well-structured representations are easily intelligible, im- 
prove information processing and can be universally shared among different NoC 
design tools. 



5 The IP Mapping Problem 

The platform-based design methodology for SoC encourages the reuse of compo- 
nents to increase reusability and to reduce the time-to-market of new designs. The 
designer of NoC-based systems faces two main problems: selecting the adequate 
set of IPs that optimize the execution of a given application and finding the best 
physical mapping of these IPs into the NoC structure. 

The main objective of the IP assignment stage is to select, from the IP repository, 
a set of IPs that minimize the NoC consumption of power, area occupied and execu- 
tion time. At this stage, no information about physical allocation of IPs is available 
so optimization must be done based on TG and IP information only. So, the result 
of this step is the set of IPs that maximizes the NoC performance. The TG is then 
annotated and an APG is produced, wherein each node has an IP associated with it. 

Given an application, described by its APG, the problem that we are concerned 
with now is to determine how to topologically map the selected IPs onto the net- 
work, such that the objectives of interest are optimized. Some of these objectives 
are: latency requirements, power consumption of communication, total area occu- 
pied and thermal behavior. At this stage, a more accurate execution time can be 
calculated taking into account of the distance between resources and the number of 
switches and links crossed by a data package along a path. The result of this pro- 
cess should be an optimal allocation of the one of the prescribed IP assignments, 
selected in an earlier stage, to execute the application, described by the TG, on the 
NoC structure. 



6 The Choice of Optimization Objectives 

Different objectives may be considered in the IP mapping problem. If the improve- 
ment of an objective leads to the deterioration of another (e.g. maximizing clock 
frequency increases power consumption), the objectives are said to be concurrent. 
On the other hand, if the improvement of an objective leads to the improvement of 
another, the objectives are said to be collaborative. Optimization problems with con- 
current and collaborative objectives are called Multi-objective Optimization Prob- 
lems (MOPs). In such problems, all collaborative objectives should be grouped and 
a single objective among those should be used in the optimization process, which 
achieves also the optimization of all the collaborative objectives in the group. How- 
ever, concurrent objectives need all to be considered in the optimization process. 
The best solution for a MOP is the solution with the adequate trade-off between 
concurrent and collaborative objectives. 

Table Q]helps choosing the minimal set of objectives to be considered in IP map- 
ping optimization stage. A up/down arrow in entry for objectives i x j means that 
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Table 1 Concurrent and Collaborative Objectives 
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an increase/reduction with respect to objective i also leads to and increase/reduction 
with respect to objective j. 

For instance, the last column of Table Q] which characterizes objective j^PE 
(i.e. number of processor elements), indicates that a reduction with respect to this 
objective yields a reduction with respect to both area and cost. Therefore, those 
three objectives are considered collaborative. This is the same case for objective 
time and clock frequency . However, the penultimate column of Table [TJ that char- 
acterizes objective power, indicates that a reduction in power leads to an increase 
in both time and clock frequency. Note that objective power must be minimized. 
Therefore, objective power is considered concurrent with both objective time and 
clock frequency . As a conclusion, the adequate trade-off can be achieved using only 
minimization functions of objectives area, execution time and power consumption. 



7 Multi-objective Evolutionary Approach 



The search space for a "good" IP mapping for a given application is defined by 
the possible combinations of IP/tile available in the NoC structure. Assuming that 
the mesh-based NoC structure has N x N titles and there are at most Af 2 IPs to 
map, we have a domain size of A^ 2 !. Among the huge number of solutions, it is 
possible to find many equally good solutions. In huge non-continuous search space, 
deterministic approaches do not deal very well with MOPs. The domination concept 
introduced by Pareto lfT6l is necessary to classify solutions. In order to deal with 
such a big search space and trade-offs offered by different solutions in a reasonable 
time, a multi-objective evolutionary approach is adopted. 

The core of the proposed tool offers the utilization of two well-known and well- 
tested MOEAs: NSGA-II Q and microGA [1]. Both adopt the domination concept 
with a ranking schema for classification. The ranking process separates solutions 
in Pareto fronts where each front corresponds to a given rank. Solutions from rank 
one, which is the Pareto-optimal front) are equally good and better than any other 
solution from Pareto fronts of higher ranks. 
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NSGA-II features a fast and elitist ranking process that minimizes computational 
complexity and provides a good spread of solutions. The elitist process consists 
in joining parents and offspring populations and diversity is achieved using the 
crowded-comparison operator |2||. MicroGA works with a very small population 
(3 to 5 individuals), which makes it very fast. A bigger population is stored on a 
population memory divided in replaceable and non-replaceable areas. The elitist 
process consists of storing the best solutions on a external memory and diversity is 
achieved using an adaptive grid Q). 

The basic workflow of both algorithms is almost the same. They start with a 
random population of individuals, where each individual represents a solution. Each 
individual is associated with a rank. The selection operator is applied to select the 
parents. The parents pass through crossover and mutation operators to generate an 
offspring. A new population is created and the process is repeated until the stop 
criterion is satisfied. 



7.1 Representation and Genetic Operators 

The individual representation is shown in Fig. [4] The tile indicates information on 
the physical location on which a gene is mapped. On a N x N regular mesh, the 
tiles are numbered successively from top-left to bottom-right, row by row. The row 
of the i th tile is given by \i/N~\ , and the corresponding column by i mod N. 



Gene 



Gene Gene ' 1 Gene 

■ ■ ■ i 

1 2 n-1 




^ "\ 






Node 
Type 


IP 
ID 


Proc. 
ID 


Tile 





Fig. 4 Chromosome encoding of an IP mapping 



The crossover and mutation operators were adapted to the fact that the set of 
selected IPs can not be changed as we have to adhere to the set of prescribed IP 
assignments. For this purpose, we propose a crossover operator that acts like a shift 
register, shifting around a random crossover point and so generating a new solution, 
but with the same set of IPs. This behavior does not contrast with the biological 
inspiration of evolutionary algorithms, observing that certain species can reproduce 
through parthenogenesis, a process in which only one individual is necessary to 
generate an offspring |fT9ll . 

The mutation operator performs an inner swap mutation, where each gene re- 
ceives a random mutation probability, which is compared against the system muta- 
tion probability. The genes with mutation probability higher than the system's are 
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Fig. 5 Application of the proposed shift crossover 
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Swap operation 

Fig. 6 Application of the proposed inner swap mutation 

swapped with another random gene of the same individual, instead of selecting a 
random IP from the repository. This way, it is possible to explore the allocation 
space preserving any optimization done in the IP assignment stage. The crossover 
and mutation strategies adopted in the IP mapping stage are represented in Fig. \5\ 
and Fig. |6l respectively. 

8 Objective Function 

During the evolutionary process, the fitness of the individuals with respect to each 
one of the selected objectives (i.e. area, time, and power) must be efficiently com- 
puted. 



8.1 Area 

In order to compute the area required by a given mapping of the application in 
question, it is necessary to know the area needed for the selected processors and 
that ocupied by the used links and switches. As a processor can be responsible for 
more than one task, each APG node must be visited in order to check the processor 
identification in the appropriate XML element. Grouping the nodes with the same 
processorlD attribute allows us to implement this verification. The total number of 
links and switches can be obtained through the consideration of all communication 
paths between exploited tiles. Note that a given IP mapping may not use all the 
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available tiles, links and switches that are available in the NoC structure. Also, ob- 
serve that a portion of a path may be re-used in several communication paths. 

In this work, we adopted a fixed route strategy wherein data emanating from 
tile i is sent first horizontally to the left or right side of the corresponding switch, 
depending on the target tile position, say j, with respect to i in the NoC mesh, until 
it reaches the column of tile j, then, it is sent up or down, also depending on the 
position of tile j with respect to tile i until it reaches the row of the target tile. Each 
communication path between tiles is stored in the routing table. The number of links 
in the aforementioned route can be computed as described in EquationQ] This is also 
represents the distance between tiles i and j and called the Manhattan distance IfTTII . 

nLinks(i,j) = \ \i/N~\ - \j/N] \ + \i mod N -j mod N\ (1) 

In the purpose of computing efficiently the area required by all used links and 
switches, an APG can be associated with a so-called routing table whose entries 
describe appropriately the links and switches necessary to reach a tile from another. 
The number of hops between tiles along a given path leads to the number of links be- 
tween those tiles, and incrementing that number by 1 yields the number of traversed 
switches. The area is computed summing up the areas required by the implementa- 
tion of all distinct processors, switches and links. 

Equation [2] describes the computation involved to obtain the total area for the 
implementation a given IP mapping M, wherein function Proc(.) provides the set 
of distinct processors used in APGm and area p is the required area for processor/?, 
function Links(.) gives the number of distinct links used in APGm, M is the area 
of any given link and A s is the area of any given switch. 

Area(M) = V] area p + (Ai + A 8 ) x Links(APG M ) + A s (2) 

peProc{APG M ) 



8.2 Execution Time 

To compute the execution time of a given mapping, we consider the execution time 
of each task of the critical path, their schedule and the additional time due to data 
transportation through links and switches along the communication path. The crit- 
ical path can be found visiting all nodes of all possible paths in the task graph and 
recording the largest execution time of the so-called critical path. The execution 
time of each task is defined by the taskTime attribute in TG. Links and switches can 
be counted using the routing table. We identified three situations that can degrade 
the implementation performance, increasing the execution time of the application: 

1. Parallel tasks mapped into the same tile: A TG can be viewed as a sequence 
of horizontal levels, wherein tasks to the same level may be executed in paral- 
lel (at the same time) allowing for a reduction of the overall execution time of 
the application. For instance, Figure [8] shows a TG that can be viewed as a se- 
quence of 7 levels: level® = {to}; leveli = {£i, £2}; level?, = {£3, tA);level-z = 
{t 5 ,t 6 ,t 7 ,t s }; leveU = {*9,*io,iii,ti2}; level 5 = {£13, £14}; and level 6 = 
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{£15}. When parallel tasks are assigned in the same processor, which also means 
that these occupy the same tile of the NoC, they cannot be executed in parallel. 

2. Parallel tasks with partially shared communication path: When a task in a tile 
(source) must send data to supposedly parallel tasks in different tiles (targets) 
through the same initial link, data to both tiles cannot be sent at the same time. 
For instance, using the mesh-base NoC of Figure Q] if the task on the most left 
upper tile have to send data to its right neighbor tile and to that on the most right 
upper tile at the same time, the initial link is common to both communication 
paths and so no parallelism can occur in such a case. 

3. Parallel tasks with common target using the same communication path: When 
several tasks need to send data to a common target task, one or more shared links 
along the partially shared path would be needed simultaneously, the data from 
both tasks must then be pipelined and so will not arrive at the same time to the 
target task. For instance, using the mesh-base NoC if FigureQ] if the task on the 
most left upper tile and that on the most right upper tile have to send data to the 
center tile at the same time, they would send it to the right and left, respectively, 
and the upper center switch would buffer the data and send it in a pipelined 
manner to the center tile. 

Equation [3] is computed using a recursive function that implements a depth-first 
search, wherein function Paths(.) provides all possible paths of a given APG and 
timet is the required time for task t. After finding the including the total execution 
time of the tasks that are traversed by the critical path, the time of parallel tasks 
executed in the same processor need to be accumulated too. This is done by func- 
tion SameProcSameLevel(.). The delay due to data pipelining for tasks on the same 
level is added by SameSourceCommonPath(.). Last but not least, the delay due to 
pipelining data that are emanating at the same time from several distinct tasks yet 
for the same target task is accounted for by function DiffSrcSameTgt{.). 

Time(M) = max > timet 

r£Paths(APG M ) \*-* t 
\t£r 

+SameProcSameLevel(r, TG, APGm) 
+SameSrcCommonPath(r, TG, APGm) 
+DiffSrcSameTgt(r, TG, APG M )) 

(3) 

Function SameProcSameLevel{.) compares tasks of a given same level that are 
implemented by the same processor and returns the additional delay introduced in 
the execution of those tasks. AlgorithmQ]shows how function SameProcLevel(.), 
that uses information from path r, application task graph T and its corresponding 
characterization graph C to compute the delay in question. 

Function SameSourceCommonPath(.) computes the additional time due to par- 
allel tasks that have data dependencies on tasks mapped in the same source tile 
and yet these share a common initial link in the communication path. Algorithm[2] 
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Algorithm 1. SameProcSameLevel(r, T, C) 

1. time :— 

2. for all t G r do 

3. forallnGTdo 

4. If T.level(t) = T.level(n) then 

5. if C. process or (t) = C .processor (n) then 

6. time := time + n.taskTime 

7. return time 



shows the details of the delay computation using information from path r, applica- 
tion task graph T and its corresponding characterization graph C. In that algorithm 
T.targets(t) yields the list of all possible target tasks of task t, A.initPath(src, tgt) 
returns the initial link of the communication path between tiles src and tgt and 
penalty represents a time duration needed to data to cross safely from one switch to 
one of its neighbors. This penalty is added every time the initial link is shared. 



Algorithm 2. SameSrcCommonPath(r, T, C) 



1. penalty :— 

2. for all t e r do 

3. iST.targets(t) > 1 then 

4. for all n G T.targets{t) do 

5. for all ri G T.targetsit) \ n / n do 

6. if C.initPath(t,n) = C.initPath(t,n') then 

7. penalty :— penalty + 1 

8. return penalty 



Function DiffSrcSameTgt(.) computes the additional time due to the fact that par- 
allel tasks producing data for the same target task need to use simultaneously at least 
a common link along the communication path. Algorithm[3] shows the details of the 
delay computation using information from path r, application task graph T and its 
corresponding characterization graph C. In that algorithm, CPath(src, tgt) is an 
ordered list w of all links crossed through src task to tgt task and penalty has the 
same semantic as in the Algorithm[2] 

8.3 Power Consumption 

The total power consumption of an application NoC-based implementation consists 
of the power consumption of the processors while processing the computation per- 
formed by each IP and that due to the data transportation between the tiles. The for- 
mer can be computed summing up attribute taskPower of all nodes of the APG and 
the latter is the power consumption due to communication between the application 
tasks through links and switches. The power consumption due to the computational 
activity is simply obtained summing up atribute taskPower of all nodes in the APG 
and is as described in Equation [4] 
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Algorithm 3. DiffSrcSameTgt(r, T, C) 



1. 


penalty :— 


2. 


for all t 6 r do 


3. 


for all t' £r \t' ^tdo 


4. 


if T.level(t) = T.level(t') then 


5. 


for all n S T.targets(t) do 


6. 


for all n' G T.targets(t') do 


7. 


if n = n' then 


8. 


w :- C.Path(t,n) 


9. 


w' :=C.Path(t',n') 


0. 


for i = to min(w. length, w' .length) do 


1. 


if to(i) = ty'(j) then 


2. 


penalty := penalty + 1 



13. return penalt y 



Power 'p(M) = > power \ (4) 

The power consumption due to communication is a very important factor and 
must be considered in order to achieve low power consumption NoC designs. An 
energy model for one bit consumption is used to compute the total energy consump- 
tion for the whole communication involved during the execution of an application 
on the NoC platform. The bit energy (Eut)> energy consumed when a data of one 
bit is transported from one tile to any of its neighboring tiles, can be obtained as in 
Equation 21 

Eut = E Sblt + E Lbit (5) 

wherein Es ut and Er Jbit represent the energy consumed by the switch and link tying 
the two neighboring tiles, respectively Q. 

The total power consumption of sending one bit of data from tile i to tile j can 
be calculated considering the number of switches and links the bit passes through 
on its way along the path, as shown in Equation[6] 

E l b f t = nLinks(i,j) x E Lbit + (nLinks(i,j) + 1) x E Sbit (6) 

wherein function nLinks(.) provides the number of traversed links (and switches 
too) considering the routing strategy used in this work and described earlier in this 
section. The function is is defined in EquationQ] 

Recall that the application TG gives the communication volume (V(t, £')) in 
terms of number of bits sent from the task t to task t' passing through a direct 
arc d t .t'- Assuming that the tasks t and t' have been mapped onto tiles i and j re- 
spectively, the communication volume of bits between tiles i and j is then V(i,j) 
= V(dt,t')- The communication between tiles i and j may consist of a single link 
lij or by a sequence of to > 1 links li. Xo , l Xo ,xn fci,x 2 ) ■ • ■ j ^ m _i,j- F° r instance, 
to establish a communication between tiles and 8 (upper most left most and lower 
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most right most) in the NoC structure of FigureQ] the sequence of links is therefore 

'0,1, h,2> ^2,5 an d 15,8- 

The total network communication power consumption for a given mapping M is 
given in Equation |7J wherein Tar gets t provides all tasks that have a direct depen- 
dency on data resulted from task t and Tile t yields the tile number into which task 
t is mapped. 



Power C {M) 



t e APG M , 

Vt' 6 Targetst 



Tr/ i \ j-iTilet.Tile+i 



(7) 



Recall that the total power consumption of the application NoC-based implemen- 
tation is given by the sum of the power consumption due the computational side of 
the application and that due to the communication involved between tiles, as explic- 
itly shown in Equation [8] 



Power (M) = Power P {M) + Power C {M) 



(8) 



9 Results 

First of all, to validate the implementation of both algorithms, these were submit- 
ted to solve mathematical MOPs and the results were compared with their original 
results. Fig. [7] shows results of both algorithms for a two objective optimization 
problem called KUR, proposed by Kursawe and used by Deb and Coello to validate 
NSGA-II @ and microGA HI, respectively. 
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Fig. 7 Results for KUR function 



Both algorithms converged to the true Pareto-front. As expected, NSGA-II found 
a higher diversity of solutions while microGA was much more faster. The parame- 
ters used for these tests were the same used by their original authors. 
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For NoC optimization, only the individual representation and the objective func- 
tions were changed, keeping the ranking, selection, crossover and mutation oper- 
ators unchanged. Different TGs generated with TGFF [4] and from E3S, with se- 
quential and parallel tasks, were used. 

Many simulations were performed to find out the setting up of the parameters 
used in NSGA-II and micro-GA. For the former, we used a population size of 600, 
mutation probability of 0.01, crossover probability of 0.8 and tournament size of 50 
and run the algorithm of 100 generations. For the latter, we used population mem- 
ory size of 1000, replaceable fraction of 0.7, non-replaceable fraction of 0.3, micro 
population of 4 individuals, mutation probability of 0.02, crossover probability of 
0.09, tournament size of 2, external memory of 200, nominal convergence of 4, re- 
placement cycle of 100, bisection of 5, and run the algorithm for 3000 generations. 

The application, represented as a TG in Fig. [8j was generated with TGFF H. 
This TG presents five levels of parallelism, formed by the mirrored left-right side 
nodes. 

Analyzing the results obtained from the first simulations, we observed that in 
order to achieve the best trade-off, the system allocated many tasks for the same 




Fig. 8 Task graph with five levels of parallelism 
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processor, which reduces area and execution time but generates hot spots iTZTll . A hot 
spot is an area of high activity within a silicon chip. Hot spots can damage a silicon 
chip and increases power consumption because of Avalanche Effect; an effect that 
increases reverse current in semiconductor materials because of temperature rising 
and other factors Ifl2l . In order to avoid the formation of hot spots, a maximum tasks 
per processor constraint was imposed in the evolutionary process. This parameter 
is decided by the NoC designer based on some extra physical characteristics. We 
adopted a maximum of 2 tasks per processor. 

The IP assignment ifTTl of the TG represented in Fig. [8] was able to discover 97 
distinct optimal IP assignments. From those 97 distinct of IP assignments, 142 op- 
timal mappings were generated. The 15 most significant solutions from the Pareto- 
front, with respect to each of the considered objectives, are listed in Table [2] and 
Table[3] Table[2]presents the IPs assigned to each node of the TG and the respective 
fitness in terms of to each of the selected objectives after running the IP assignment 
step. Table[3]presents the tile where each assigned IP was mapped and Table|4]shows 
the respective fitness in terms of each of the selected objective after completing the 
IP mapping step. The first five solutions impose shorter execution times; the next 
five require smaller hardware areas and the last five solutions present lower power 
consumptions. The differences in execution time and power, observed when com- 
paring the data from both tables is because of the inflicted penalty in execution time 
and power due to the use of shared links and switches. 

Fig. |9]-(a) represents the time x area trade-off, Fig. |9]-(b) depicts the power x 
time trade-off and Fig.|9]-(c) plots the power x area trade-off. As we can observe, 
comparing the dots against the line of interpolation, the trade-off between time and 
area and between power and time is not so linear as the trade-off between power and 
area. Fig.|9]-(a) shows that solutions that require more area tend to spend less execu- 
tion time because of the better distribution of the tasks allowing for more parallelism 
to occur. Fig. |9]-(b) shows that solutions that spend less time of execution tend to 
consume more power because of IP's features, such as higher clock frequency, and 
physical effects like intensive inner-electrons activity. Fig. |9]-(c) shows a linear re- 
lation between power consumption and area. Those values and units are based on 
E3S Benchmark Suite 0. 

Figure [TOl-(a) shows the Pareto-front discrete points. Figure [TOl-(b) shows the 
Pareto-front formed by the Pareto-optimal solutions. Note that many solutions have 
very close objectives values. 

For a TG of 16 tasks, a 4 x 4 mesh-based NoC is the maximal physical structure 
necessary to accommodate the corresponding application. Table [4] shows that no 
solution used more than ten resources to map all tasks. The unused 6 tiles may 
denote a waste of hardware resources, which consequently lead to the conclusion 
that either the geometry of the NoC is not suitable for this application or the mesh- 
based NoC is not the ideal topology for its implementation. 
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Fig. 9 Trade-offs representation of the 142 optimal IP mappings obtained for the task graph 
of Fig. ! 
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Table 4 Characteristics of optimal IP mappings for the task graph of Fig. [8] 



solution 


time (s) 


area 


power (W) 


Used tiles 


1 


0.1225 


19.9514 


52.9501 


10 


2 


0.129 


6.386 


21.1251 


9 


3 


0.1362 


6.506 


21.1251 


9 


4 


0.144 


6.2951 


23.5501 


8 


5 


0.1503 


6.566 


19.6001 


9 


6 


0.2073 


3.6056 


15.8401 


8 


7 


0.2583 


3.787 


15.0401 


8 


8 


0.1503 


4.4066 


18.7001 


8 


9 


0.1829 


4.5676 


16.8701 


9 


10 


0.1819 


4.6276 


17.2701 


9 


11 


0.2073 


4.267 


15.0401 


8 


12 


0.2073 


3.7856 


15.8401 


8 


13 


0.1819 


4.8076 


16.8701 


9 


14 


0.1819 


4.6276 


17.2701 


9 


15 


0.1503 


4.6466 


18.7001 


8 



1 (xl0" 5 m 2 ). 



Table 5 Processors of solution #8 from Tableland Table|3] 



TG Node 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


Proc ID 


32 


32 


15 


13 


17 





6 


17 


30 


6 


13 





30 


15 


23 


23 


IP ID 


942 


937 


458 


378 


490 


43 


240 


480 


855 


216 


379 


13 


862 


456 


724 


719 


Tile 








4 


5 


10 


6 


1 


10 


9 


1 


5 


6 


9 


4 


8 


8 



As a specific mapping example, we detail solution #8 in Table[3] which seems to 
be a moderate solution with respect to every considered objectives. Table|5]specifies 
the processors used in solution #8 and Fig. QT| shows the mapping of this solution 
into the mesh-based NoC. We can observe that all parallel tasks were allocated in 
the distinct processors, which reduces execution time. The number of processors 
were minimized based on the optimization of the objectives of interest and this min- 
imization was controlled by the maximum tasks per processor constraint to avoid 
hot spots |2T1 . The processors were allocated in such way to avoid delay of commu- 
nication due to links and switches disputed by more than one resource at the same 
time. 
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Fig. 11 Mapping of solution #8 from Tableland Table|3] 



10 Conclusions 



The problem of mapping IPs into a NoC topology is NP-hard J6) key research prob- 
lems in NoC design 1151 . In this paper we propose a decision support system based 
on MOEAs to help NoC designers allocate a prescribed set of IPs into a NoC phys- 
ical structure. The use of two different MOEAs consolidates the obtained results. 

Structured and intelligible representations of a NoC, a TG and of a repository 
of IPs were used and these can be easily extended to different NoC applications. 
Despite of the fact that we have adopted E3S Benchmark Suite as our repository 
of IPs, any other repository could be used and modeled using XML, making this 
tool compatible with different repositories. 

The proposed shift crossover and inner swap mutation genetic operators can be 
used in any optimization problem where no lost of data from a individual is ac- 
cepted. 

Future work is three-fold: adopting a dynamic topology strategy to try to evolve 
the most adequate topology for a given application; exploring the use of differ- 
ent objectives based on different repositories and proposing an interfacing mecha- 
nism with a hardware description simulator to integrate our tool to the NoC design 
platform. 
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1 Introduction 

In our society, various combinatorial optimization problems exist and we 
must often solve them, for e.g. scheduling, delivery planning, circuit design, 
and computer wiring. Then, one of the important issues in science and engi- 
neering is how to develop effective algorithms for solving these combinatorial 
problems. 

To develop effective algorithms for such combinatorial optimization prob- 
lems in the real world, the standard combinatorial optimization problems are 
intensively studied: for example, traveling salesman problems, quadratic as- 
signment problems, vehicle routing problems, packet routing problems, and 
motif extraction problems. Among them, the traveling salesman problem 
(TSP) is one of the most standard combinatorial optimization problems . 
The TSP is described as follows: when a set of N cities and distances dij 
between cities i and j are given, find an optimal solution, or the shortest- 
length tour visiting all the cities once. Namely, the goal of TSP is to find a 
permutation cr of the cities that minimizes the following objective function: 

iV-l 
L{<t) = 2_^ d a {k)a{k+l) + d a (N)<j(l), (1) 

fc=l 
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where L(cr) is the tour length of the TSP with cr. If dy = dj% for all i and j, 
the TSP is symmetric; otherwise, it is asymmetric. In this chapter, we deal 
with the symmetric TSP. 

For an TV-city symmetric TSP, the number of all possible tours is (N — 
l)!/2. Thus, the number of tours factorially diverges if the number of cities 
increases. It is widely acknowledged that the TSP belongs to a class of NP- 
hard. Therefore, it is required to develop an effective approximate algorithm 
for finding near-optimal solutions in a reasonable time frame. 

To discover approximate solutions, various heuristic algorithms have al- 
ready been proposed. In 1985, Hopficld and Tank proposed an approach for 
solving TSPs by using a recurrent neural network. They applied descent 
downhill dynamics of the recurrent neural network to obtain approximation 
solutions of TSPs |T . Although this approach is interesting from the view- 
point of an application of neural dynamics to real engineering problems, such 
as combinatorial optimization, it has two drawbacks. 

The first drawback is that this approach has a local minimum problem: it 
uses simple descent downhill dynamics of the neural network to obtain better 
solutions of TSPs; states of the neural networks can be stuck at undesirable 
local minima. To resolve the local minimum problems, two main strategies 
that uses chaotic dynamics have been proposed. The first solution is to inject 
chaotic noise into the dynamics of the neural network [2,3,4,5,6 . The second 
solution is to replace the descend downhill dynamics with chaotic dynamics 
[TUHHHlIin]- Recently, these strategies have been applied to combinatorial 
optimization in engineering, for example, channel assignment problems |llj . 
frequency assignment problems [12 , multicast routing problems |13j . and 
broadcast routing problems [13j . The performance of using chaotic dynamics 
shows that the algorithm finds an optimal or near-optimal solutions of the 
problems. In this chapter, we review basic theories of these strategies in 
Sects. O and O 

Although the first drawback can be resolved by the above-mentioned 
strategies, namely chaotic noise injection and chaotic dynamics in recurrent 
neural networks, there still exists the second drawback. The methods based 
on the recurrent neural networks cannot be applied to large scale instances 
of combinatorial optimization, because it takes a huge amount of memories 
to construct the neural network. In addition, it is not so easy to obtain a 
feasible solution. In the method, a solution of TSPs is encoded by a firing 
pattern of the neural network. Thus, a solution is generated only in the case 
that the firing pattern of the neural network completely satisfies a constraint. 

To resolve the second drawback, an approach in which heuristic algorithms 
are controlled by the chaotic dynamics has been proposed [T3HT21III3III1IIE] ■ 
In this approach with chaotic dynamics, to avoid local minima, execution 
of a heuristic algorithm, such as the 2-opt algorithm , the 3-opt algorithm , 
the Or-opt algorithm [!§] , the Lin-Kernighan (LK) algorithm 20J , and the 
stem-and-cycle (S&C) ejection chain method [2"T1|2"2"] . is controlled by chaotic 
dynamics. 
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In [2J[T31[Tni[Hl[TH], to generate the chaotic dynamics, a chaotic neural 
network [231211 is used. In the chaotic neural network, the basic element is 
a chaotic neuron proposed by Aihara et al. (2I3I23- The model introduced 
two important properties which real nerve cells have, namely graded response 
and refractoriness . The latter means that when a neuron has just fired, the 
firing of this neuron is inhibited for a while by the refractoriness. 

In [TJJ[TS1[TH1[T71[IB] , execution of the local search algorithm is encoded by 
firing of the chaotic neuron. If the chaotic neuron fires, the corresponding 
local search algorithm is executed. Because the firing of the chaotic neuron 
is inhibited by the refractoriness, frequent firing of the chaotic neuron, or 
frequent execution of the local search algorithm is restricted. This can be a 
mechanism for the chaotic search to escape from local minima efficiently. It 
is reported that the refractoriness implemented in the chaotic neuron model 
leads to a higher solving ability than the tabu search which has a similar 
strategy of searching solutions as the chaotic search [15] . 

On the basis of the above idea that the refractoriness of the chaotic neuron 
can be used for an effective escape from local minima, a chaotic search method 
which controls the 2-opt algorithm has already been proposed P31|T5]. Al- 
though the 2-opt algorithm is the simplest local search algorithm, the chaotic 
search method with the 2-opt algorithm shows good results. In [14], in the case 
of solving an A^-city TSP, N x N chaotic neurons are prepared and arranged 
on N x N grid. Here, N x N neurons correspond to the number of possible 
ways for constructing a new tour by the 2-opt improvement. As a result, this 
chaotic search method shows solving performance with less than 0.2% gaps 
from the optimal solution for instances with the order of 10 2 cities [14] . 

On the other hand, the tabu search [2"5lT2l)] is also a quite effective strategy 
for escaping from local minima. The chaotic search method in }15| is based 
on the tabu search which is implemented on a neural network, because both 
the tabu search and the chaotic search have similar strategies that forbid 
backward moves. In [15] . by assigning one neuron to one city, only N chaotic 
neurons are used to solve an JV-city TSP. As a result, this chaotic search 
method can be applied to large scale examples, such as the 85,900-city one, 
and it has high solving performance with less than 3.0% gaps from the optimal 
solutions for instances with the order of 10 4 cities [15 . In addition, this 
strategy has been improved by introducing other local searches, such as the 
Or-opt algorithm [19] , the Lin-Kernighan (LK) algorithm [20] , and the stem- 
and-cycle (S&C) ejection chain method [2TH22]. In [HHSHUMlISniEIlIITlEl] , 
these sophisticated local searches are controlled by chaotic dynamics. The 
results show that large scale TSPs could be solved by these chaotic searches. 

In this chapter, the chaotic search that resolves the second drawback in the 
Hopfield-Tank neural network approach is described in In [3l we also review 
several applications of this chaotic search to real engineering problems: the 
vehicle routing problems in l3.U the motif extraction problems in !3.2l and the 
packet routing problems in[ 



134 T. Ikeguchi et al. 

2 Methods 

In this section, we review the methods for solving combinatorial optimization 
problems by chaotic dynamics. 

2.1 Chaotic Noise 

Performance improvement by the chaotic noise 

Effectiveness of the chaotic noise for shaking the states of the solution search 
has been shown by many papers 2,3,4,5,6. In Figs. Q] and[21 the solvable per- 
formances of the Hopficld-Tank neural networks pQ with stochastic noise and 
chaotic noise are compared by applying to the TSP and QAP, respectively. 
For these comparisons, the logistic map chaos, z(t + 1) = az(t)(l — z(t)), is 
introduced as a most simple example of chaotic sequence, with the parame- 
ters a — 3.82, a — 3.92, and a — 3.95. As the stochastic noise, the Gaussian 
white noise is introduced. The horizontal axis is noise amplitude (3 and the 
vertical axis is the percentage of the optimum solution obtained by 1,000 dif- 
ferent initial conditions. From the figures, it is clear that the chaotic noise is 
effective for combinatorial optimization algorithm using the recurrent neural 
networks. Although the original Hopfield-Tank neural network quickly con- 
verges to a stable state corresponding to a local optimal and does not offer 
the optimum solution for these problems in any cases, its performance can be 
much improved by adding the chaotic noise, which makes the solvable perfor- 
mances almost 100% for the 20-city TSP as shown in Fig. [T] and around 20% 
for the 12-node QAP as shown in Fig. respectively. The white Gaussian 
noise also improves the performance of the Hopfield-Tank neural network, 
but its solvable performances are around 80% and 6% for the TSP and the 
QAP, respectively, and much lower than the chaotic noise cases. These results 
show that the chaotic noise has much better solvable performance than the 
stochastic noise. 

The Hopfield-Tank neural network can be applied to combinatorial opti- 
mization problems, based on its minimization property of the energy function 

-. n n n n 

i=l j=l k=l 1=1 i=l k=l 

which always decreases by asynchronous update of each neuron by the fol- 
lowing simple equation: 



X ik (t + l) = f 






j=l 1=1 



(3) 



where Xit(t) is the output of the ifcth neuron at time t, Wikji is the connection 
weight between the ikih and j'/th neurons, and dik is the threshold of the ifcth 
neuron. 
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Fig. 1 Solvable performance of the recurrent neural networks with the chaotic noise 
and the white Gaussian noise on a 20-city TSP 
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As already been applied in many previous researches, the energy function 
for solving the TSPs [Iim i8l[l"0ll2ll3] can be defined by the following equation: 



Etsv — A 



N / N 



E 5>**(*)-i 



N N N 




N ( N \ ' 



+ B EEE ^^ife(*){a;jfc+i(*) + ^fc-i(t)}, 

i=i j=i fe=i 



(4) 



where Af is the number of cities, dy is the distance between the cities i and 
j, and A and -B are the weight of the constraint term (formation of a closed 
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tour) and that for the objective term (minimization of total tour length), 
respectively. In this neural network, firing of the ikth neuron means that the 
city i is visited at the k order. From fl2} and (JU, the connection weights Wijki 
and the threshold Oijki can be obtained as follows: 

Wikji = -A{Sij(l - 5 U ) + S k i(l - Sij)} - Bdij(8i k+1 + 8i-k-i), (5) 

% - 2A, (6) 

where Sij = 1 if i = j, Sij = otherwise. Using these connection weights and 
the thresholds for the update equation in (|3|), the better solution of the TSP 
may be found by simple neuronal updating, that appears as a firing pattern. 
For the QAPs whose objective function is 

N N 

F (p) = EE^ 6 pWp(j)> (7) 

which has to be minimized by finding the optimal permutation p, we use the 
following energy function: 



Eqap = A 



N / N \ I \ N / N \ • 



N N N N 

4=1 j=l fc=l 1 = 1 

where N is the size of the problem, ay and bki are N x N matrices given in 
the problem definition in (J7J), and A and B are the weight of the constraint 
term (making p a permutation) and that for the objective term (minimization 
of the objective function), respectively. In this neural network formulation, 
firing of the ifcth neuron means that the element i is assigned to the fcth 
location of the permutation. By transforming (JSJ) to the form of the energy 
function of the recurrent neural network in |2j , the connection weights and 
the thresholds for solving the QAPs are obtained as follows: 

mkji = -A{5ij(l - Su) + S k i(l - S^)} - Bdijbki, (9) 

% = 2A. (10) 

Using these connection weights and the thresholds for neuronal updates by 
(J3J), the better solutions of the QAP may be obtained. 

However, the Hopfield-Tank neural network is well known to have the lo- 
cal minimum problem, because the energy function simply decreases and it 
can search simply one solution in a huge number of local optimal solutions. 
To improve the performance of such a neural network, the chaotic noise or 
stochastic noise has been applied to this optimization neural network to im- 
prove the solutions by avoiding trapping at the undesirable optimal states. 
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In [3H2HS], the following update equation is used to apply the noise to the 
neural network, 



Xik(t+ 1) = / 



N N 

y^ y^ Wikjixji(t) + Oik + f3z ik {t) 



(ii) 



where Zik(t) is a noise sequence added to the ikth neuron, f3 is the amplitude 
of noise, and / is the sigmoidal output function, f(y) = 1/(1 + cxp(— y/e)), 
respectively. The noise sequence introduced as Zik{t) should be normalized 
to zero mean and unit variance. 

The results in Figs. Q] and [5] are obtained by using the neural network with 
noise described above, with changing the noise amplitude /3. From the results, 
the noise amplitude /3 for the best performance differs among the noise. By 
comparing the best solvable performances in each noise, effectiveness of the 
chaotic noise can be clearly shown by such simple experiments. 

Analysis on the effects of the chaotic noise 

In order to know why the chaotic noise is more effective than other noise, 
such as the Gaussian white noise, effectiveness of the chaotic noise has been 
analyzed from various aspects. 

In [5] , Hayakawa et al. compared the performance of the neural networks 
with the chaotic noise generated by the logistic map and those with randomly 
shuffled sequences of the chaotic noise, whose temporal structure, such as 
auto-correlation, is broken. Their results show that the performance with 
the random shuffled sequence becomes much worse than the original chaotic 
sequence. From such results, they anticipated that the temporal structure of 
the chaotic noise is important in the chaotic search. 

In [3], Hasegawa et al. presented much clearer results showing the impor- 
tance of the temporal structure of the chaotic noise, by applying the method 
of surrogate data [33J . They introduced three algorithms for generating surro- 
gate data. The first one is the random shuffle algorithm, which is almost the 
same as the method which was introduced in [2] mentioned above. It preserves 
the empirical histogram of the original data. The second one is the Fourier 
transformed surrogate algorithm, which generates stochastic data preserving 
the auto-correlation function and the power spectrum of the original data. 
Such surrogate data can be generated by applying the discrete Fourier trans- 
form to the original data for obtaining the power spectrum, randomizing the 
phase with keeping the same power spectrum, and then applying the inverse 
discrete Fourier transform to the phase randomized spectrum. The generated 
sequence will have the same power spectrum and auto-correlation function as 
the original data. The third algorithm also preserves the auto-correlation and 
power spectrum of the original data, and additionally the empirical histogram 
of the original data. 
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The results in [5J show that the neural networks with the noise sequences 
generated by the second and the third surrogate algorithms, which preserve 
the auto-correlation function and the power spectrum of the original chaotic 
sequence, have almost the same performances as those of the neural network 
with the original chaotic sequence. This result clearly shows that temporal 
structure, such as auto-correlation function, of each noise is an important 
factor for high performance of the chaotic search. 

In [3] , Hasegawa and Umeno focused on the shape of the auto-correlation 
function of the chaotic sequences which leads the neural network with high 
performance. Such chaotic sequences have the auto-correlation that gradually 
converges to zero with damped oscillation. Such chaotic sequences with sim- 
ilar auto-correlation have also been utilized in the chaotic CDMA researches 
[5U[55]. I n CDMA, minimization of the cross-correlation among the spread- 
ing sequences leads to the lower mutual interference. In the chaotic CDMA 
researches such as [34l[35], the chaotic sequences, whose auto-correlation are 
C(t) rs c x (— r) r , have been used, that is similar to the chaotic sequences 
which leads high performance on the optimization neural network described 
above, where r is the damping factor, r is the lag, and c is a constant. It 
has been shown that the auto-correlation with the r = 2 — \/3 leads to the 
minimum cross-correlation among the sequences. By using such a sequences, 
performance of the bit error rate in the CDMA communication could be 
improved in [35) . 

Hasegawa and Umeno also investigated the performance of the neural net- 
work with such noise whose auto-correlation is C{t) kcx (— r) r , and showed 
that higher performance can be realized only by the noise with positive r 
which is similar auto-correlation as the original chaotic sequence. Further- 
more, in [Bj, Minami and Hasegawa showed that injection of negative auto- 
correlation sequences leads to the lower cross-correlation that may be realized 
by the same mechanism as the chaotic CDMA [511155] . From these researches, 
it has been shown that the neural networks with the chaotic noise have higher 
solving abilities because their negative auto-correlation makes smallest cross- 
correlation between the neurons that leads the high dimensional searching 
dynamics of the neural network to the most complicated dynamics, and such 
dynamics makes the performance much higher than other noise sequences 
such as the white noise. 



2.2 Recurrent Chaotic Neural Networks 
Performance improvement by the chaotic neural networks 

The chaotic neural network [23 24.36 has an inherent property to escape from 
any fixed points except that of all the resting neurons due to accumulated 
refractory effects. This property was first applied to dynamical associative 
memory [23 , 57] , then to the optimization methods based on the Hopficld- 
Tank neural networks pQ , and effectiveness of chaotic dynamical searches has 



Theory and Applications of Chaotic Optimization Methods 139 

also been shown by many papers [TJIHlISlUni ■ The chaotic neural network is 
a neural network model consisting of the chaotic neurons, which have the 
chaotic dynamics. The chaotic neuron is an extension of the Nagumo-Sato 
binary neuron model, which has an exponentially decreasing refractory effect, 
to an analog version by replacing the Heaviside output function to the sig- 
moidal output function. The chaotic neuron map can be described as follows: 

y(t + l) = ky(t)-af(y(t))+a, (12) 

where y(t) is the internal state of the neuron at time t, k is the decay parame- 
ter of the refractory effects, a is the scaling parameter of the refractory effects, 
a is a positive bias, and / is a sigmoidal function, f(y) = 1/(1 + exp(— y/e)). 
The chaotic dynamics can be represented by this bimodal nonlinear map 

nam]- 

The chaotic neural network is a network composed of the chaotic neurons, 
which is defined as follows 23, 36l[3Z| : 



r)ij(t + 1) = kmWijit) + /]/, WijkiXki(t), (13) 

k i 
dj (t + l) = k r (ij (t) - axij (t) + dij , (14) 

Xij (t + 1) = flVij (t + 1) + Cij (t + 1)] , (15) 

where r)ij (t) is the internal state for mutual connections of the ijth neuron at 
time t, Qj(t) is the internal state for the refractory effects of the ijth neuron 
at time t, a%j is the positive bias for the ijth neuron, k m and k r are the decay 
parameters for the mutual connections and refractory effects, respectively. 

In P3HHH]; the chaotic neural network with a single internal state which 
corresponds to setting k r — k m is used. Effectiveness of such chaotic dynam- 
ics has been shown by Nozawa in [7], by comparing the performances with 
those of the stochastic models on the basis of extension of the continuous-time 
Hopficld neural network model to a discrete-time model with adding nega- 
tive self-feedback connections for each neuron. Such a negative self-feedback 
connection corresponds to the refractory effects in the above chaotic neural 
network model, and this neural network has the chaotic dynamics as well. To 
improve the performance of the chaotic search, Chen and Aihara proposed 
chaotic simulated annealing by gradually reducing the effects of the chaotic 
fluctuation in the searching dynamics [5], and showed that the performance 
can be much improved. 

Analysis on the effects of the chaotic neural network 

In the chaotic noise approach in the previous section, the chaotic dynamics 
has been introduced as additive noise to the gradient dynamics. In contrast 
with such an approach, the chaotic neural network approach uses the search- 
ing dynamics whose dynamics itself is chaotic. Therefore, the searching dy- 
namics of the chaotic neural networks has various well-known characteristics 
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of the chaotic dynamics, such as orbital instability , self-similarity , and so 
on, and has better performance than the chaotic noise approach [55]. 

Such chaotic searching dynamics of the chaotic neural networks has been 
analyzed by calculating the Lyapunov exponents [9l[T0]. Although it is not 
easy to estimate the Lyapunov exponents accurately for such high dimen- 
sional chaotic dynamical systems, clear results have been obtained. Yamada 
and Aihara calculated the Lyapunov exponents of the single internal state 
model [5] . They analyzed the performance of the chaotic neural network with 
changing its parameter values, and showed that sum of the positive Lya- 
punov exponents of the high performance chaotic dynamics becomes close to 
zero. They argued that such results mean that the edge of chaos , between 
the periodic dynamics and the chaotic dynamics, has the best performance 
for combinatorial optimization. Hasegawa et al. also analyzed the relation 
between the solvable performances of the chaotic neural networks and its 
Lyapunov exponents on the two-internal-state model described in (fT3|) and 
(fT4"|) . They showed that it is possible to tune the dynamics of the chaotic neu- 
ral network by changing the balance between two decay parameters, k r and 
k m , for the internal states rjij(t) and Qj(t), respectively, and obtained similar 
results that the chaotic dynamics with smaller positive Lyapunov exponents 
has the best performance. They also calculated the coefficient of variation, 
as a complexity index of the firing interval of the neurons, and showed that 
the chaotic dynamics with small Lyapunov exponents has higher complexity 
that makes the chaotic neural network realize high solvable performances. 

From the obtained results in various researches on the both approaches 
using chaotic dynamics, the chaotic noise and the chaotic neural networks, 
for the recurrent neural networks, described in the previous and present sec- 
tions, it has been understood that chaos makes the searching dynamics very 
complicated and the performance improved. Such chaotic dynamics is more 
complex than the stochastic dynamics in a sense, and the better performance 
can be realized. By further researches, it may be possible to completely clar- 
ify the mechanism of the effectiveness of the chaotic search in the recurrent 
neural network in the near future. 

2.3 Chaotic Search That Controls Execution of 
Heuristic Algorithm 

Recurrent chaotic neural networks are effective for solving combinatorial opti- 
mization problems as shown in the previous section l^T^l However, the method 
cannot be applied to very large instances because it needs huge amounts of 
memory to construct the neural network. In addition, it is not so easy to 
obtain feasible solutions because a firing pattern of a chaotic neural net- 
work encodes a solution. Thus, a solution is generated only in the case that 
the firing pattern of the neural network satisfies the constraints. To resolve 
these serious problems, the second approach, in which heuristic algorithms 
are driven by the chaotic dynamics, has been proposed [3]. 
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In the approach, execution of a local search algorithm is controlled by the 
chaotic dynamics. The basic element is a chaotic neuron proposed by Aihara 
et al. [53jrj3]. Execution of the local search algorithm is encoded by firing of 
the chaotic neuron. If the chaotic neuron fires, the corresponding local search 
algorithm is executed. After a neuron has just fired, the next firing of this neu- 
ron is inhibited for a while by the refractoriness of the chaotic neuron. Thus, 
frequent firing of the chaotic neuron, or frequent execution of the local search 
algorithm is restricted. Therefore, the chaotic search can escape from local 
minima efficiently. Then, it is reported that the refractoriness implemented 
in the chaotic neuron model leads to equivalent or even higher solving ability 
than tabu search which has almost the same strategy of searching solutions 
as the chaotic search. 

Using the above-mentioned idea, chaotic search methods have already been 
proposed to find the near optimal solutions or approximate solutions for 
combinatorial optimization problems such as traveling salesman problems 
[T4"ll39U40[[r5II 1 611 1 8ll2"71ITTll3"2] . quadratic assignment problems [41132"] . vehicle 
routing problems 43,44,45 , motif extraction problems 46,47,48,49,18 , and 
packet routing problems [501EI1E21ES1IM1ES1ES1EZ1ISB]- In this section, we 
describe the simplest chaotic search method for solving TSP 14,39,40,15 . 
In the method, execution of the 2-opt algorithm is driven by the chaotic 
neurodynamics . The 2-opt algorithm exchanges two paths with other two 
paths until no further improvement can be obtained (Fig. [2]). However, a tour 
obtained by the 2-opt algorithm is not a global optimum but a local optimum. 
To jump from such a local optimum, we applied the chaotic neurodynamics to 
the 2-opt algorithm 14.41:15]. To realize the chaotic search and to avoid local 
minima by the chaotic dynamics, a chaotic neuron is assigned to each city. 
Then, if a chaotic neuron fires, the local searches related to the corresponding 
city are carried out. 

Dynamics of the ith chaotic neuron is defined as follows: 



*i(t+l) = /(lfc(t+l)), (16) 

f(y) = t- — 1 - r ^rv W 

1 + exp(-2//e)) 

y l (t + i) = e l (t + i) + at + i), (is) 




Fig. 3 An example of the 2-opt algorithm. In this example, a(i) is the next city 
to i. Two paths (i-a(i) and j-a(j)) are deleted from the current tour, then new two 
paths, i-j and a(j)-a(i), are added to obtain a shorter tour 
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where Xi(t + 1) is an output of the ith chaotic neuron at time t + 1; /(•) is 
a sigmoidal function; yi(t + 1) is an internal state of the ith chaotic neuron 
at time t + 1. The internal state is decomposed into two effects: a gain effect 
and a refractory effect. 

The gain effect is defined by 

&(* + 1) - max{/3(t)zi. u (i) + Q(t)}, (19) 

j 

where Aij (t) is difference between the length of a current tour and that of a 
new tour in which city j is next of city i after applying the 2-opt algorithm 
to city i (Fig. |3J). Q(t) is a refractory effect of the jth city at time t which 
is defined in (j2"Tj) . In (|19D . the refractory effect of the jth city is included to 
avoid local minima. Let us assume that a searching state now gets stuck at 
a local minimum. Then, we calculate a maximum value of A^. In this case, 
the maximum value is Aij — because the current solution is a local optimal 
solution or an optimal solution. Thus, to select other cities, the refractory 
effect is included in (flT)l) . In this equation, j3{t) is a scaling parameter, which 
increases in proportion to time t as follows: 

0{t + 1) = /3{t) + A. (20) 

If we use (|2"0|l. the searching space is gradually limited, which has a similar 
characteristic as the simulated annealing [55]. If (3{t) takes a small value, the 
method can explore a large solution space. On the other hand, if (3(t) takes 
a large value, the method works like a greedy algorithm. 

The refractory effect works to avoid the local minima, which has a similar 
effect as a memory in the tabu search [25,26.. In the tabu search, to avoid 
a local minimum, previous states are memorized by adding them to a tabu 
list and are not allowed again for a certain temporal duration called a tabu 
tenure. In the case of the chaotic search, past firing is memorized as previous 
states to decide the present strength of the refractory effect, which increases 
just after corresponding neuron firing and recovers exponentially with time. 
Thus, while the tabu search perfectly inhibits to select the same solutions 
for the certain period, the chaotic search might permit to select the same 
solutions if a corresponding neuron fires due to a larger gain effect or an 
exponential decay of the refractory effect. The refractory effect is expressed 
as follows: 

t 
Q(t + l) = -aJ2 k rXi(t-d) + e (21) 

= krCi(t)-aXi(t) + 0(l-kr), (22) 

where a controls the strength of the refractory effect after the firing (0 < a); 
the parameter k r is the decay parameter of the refractory effect (0 < k r < 1); 
and 9 is the threshold value. If the neuron frequently fires in its past history, 
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the first term of the right hand side of (|55|) becomes negative. Then the 
neuron leads to a resting state. By using the chaotic neuron, we can realize 
the tabu search method as a special case [351 [3|)J EJUS]. The conventional 
tabu effect can be described by setting the parameters a — > oo and k r = 1 in 
the refractory effect C,[ of the chaotic neuron as follows: 

s-1 

$ = -a y £ / 4x ij (t-d) + 9, (23) 

d=0 

where s corresponds to a tabu tenure. If the neuron fired during s steps, firing 
of this neuron is inhibited by the parameters a and k r . It is considered that 
the computational cost is almost the same as the tabu search method. We 
have already shown that the chaotic search method shows better performance 
than the tabu search method in several examples [mrai^imii^MISOllHl 

152] . 

When we solve the TSP, the following procedure is used: 

1. An initial tour is constructed, for example by the nearest neighbor method. 

2. The tour is improved by the 2-opt algorithm controlled by chaotic 
dynamics. 

a. A city i is selected from the neurons whose internal state has not been 
updated yet. 

b. A city j is selected in such a way that the gain effect is maximum. 

c. If the ith neuron fires, city i and city j are connected by the 2-opt 
algorithm. 

d. The steps a)-c) are repeated until all neurons are updated. 

3. One iteration is finished. Then, the step 2 is repeated. 

To extend the above-mentioned algorithm, we have also proposed a new 
chaotic search method [17 . In the method, execution of the LK algorithm [2U] 
is controlled by the chaotic dynamics. The LK algorithm [20 is one of the 
most powerful variable depth search methods. It can explore better solutions 
than the adaptive fc-opt algorithm because the adaptive fc-opt algorithm is 
based on a simpler rule. The number of exchanged links k increases when a 
gain of a (k + l)-opt improvement is larger than that of a /c-opt improvement. 
As a result, the chaotic search method using the LK algorithm shows solving 
performance with less than 0.7% gaps from the optimal solution for instances 
with the order of 10 4 cities and can be applied to large scale instances with 
the order of 10 5 cities [17] . 

On the other hand, the S&C ejection chain method [2TII60J is also one of 
the most effective variable depth search methods. It is reported that the S&C 
ejection chain method leads to better solutions than the LK algorithm [60 . 
One of the reasons is that the S&C ejection chain method can explore more 
diversified solution space, because it introduces an S&C structure, which is 
not a tour. Namely, the S&C ejection chain method can explore infeasible 
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Table 1 The results of the local search algorithms and the chaotic search methods. 
Results are expressed by percentages of gaps between obtained solutions and the 
optimal solutions (%) 





Local search 




Chaotic search 


2-opt 


2opt + 


LK SC 


2-opt 




2opt+ 


LK SC 


Instance 


Or-opt QU 


m m\ 


hhesiidiitsi 


Or-opt |27] 


H3 m 


pcb442 7.473 


3.970 


2.298 2.148 


1.034 




0.409 


0.336 0.369 


pcbll73 9.885 


6.238 


2.903 2.680 


1.692 




0.804 


0.599 0.452 


pr2392 9.563 


5.294 


4.225 3.252 


1.952 




1.153 


0.619 0.466 


rl5915 9.395 


6.244 


4.311 3.268 


2.395 




1.291 


0.748 0.702 


rlll849 8.752 


5.567 


4.066 3.248 


2.223 




1.139 


0.708 0.652 



solution space. However, the S&C ejection chain method also gets stuck at 
local minima because it is also a greedy algorithm. 

Table [T] shows performances of chaotic search methods for benchmark in- 
stances of TSPLIB HI] (see |27017I32| for details of the algorithms). Although 
the computational cost of the chaotic search methods is larger than that of 
the local search method because the refractory effect is calculated to avoid 
local minima, the performance of the local search methods is much improved 
by the chaotic dynamics. 

3 Applications 

In the previous section, we described the three basic approaches for solving 
combinatorial optimization problems with chaotic dynamics. In this section, 
we review the application of the chaotic search methods to three important 
engineering applications: vehicle routing problems, motif extraction prob- 
lems, and packet routing problems. 

3.1 Vehicle Routing Problems 

To plan an efficient schedule of the delivery, the vehicle routing problem 
(VRP) is widely studied 62 25 ,26 ,63 64,43,44,45 . In this section, we explain 
a chaotic search method for solving VRP. Then, we show the results for 
the Solomon's benchmark problems |65) and the Gehring and Homberger 
benchmark problems |66) . 

The VRP consists of a depot, vehicles and customers. The depot is a 
departure and an arrival point of the vehicles. Each vehicle has a weight 
limit and visits the customers to satisfy their demands. The customers are 
visited only once by one vehicle. Then, the object of the VRP is to minimize 
the number of vehicles and the total travel distance. Generally speaking, a 
primary object of the VRP is to minimize of the number of vehicles. Thus, 
we use the following objective function: 
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9 (S)=^A+7Xm, (24) 

1=1 

where S is a solution (a set of tours of all vehicles), m is the number of 
vehicles, Di is the total travel distance of the Ith vehicle, and 7 is the scaling 
parameter. Because the first priority of the VRP is to reduce the number of 
vehicles, we set 7 large. We deal with a VRP with time windows (VRPTW). 
In the VRPTW, each customer has its own time window, and the vehicles 
have to visit the customers within the time window. If the time windows are 
violated, solutions are infeasible. 

In the chaotic search method, we use two simple local searches. The first 
one is to exchange the customer for another one, and the second one is to 
relocate the customer to another place. In the method of [44,43,45], 2n 
neurons are needed to solve an n-customer problem. Each neuron corresponds 
to each customer. If a neuron fires, a customer corresponding to the neuron 
is exchanged or relocated. 

To realize the chaotic search, each neuron has a gain effect and a refractory 
effect. These effects of the ijth. neuron are defined as follows: 



Z ij (t + l) = max{/3A ijs }, (25) 

s 

t 
(ij (t + l) = -aJ2 k r^j {t~d)+ 6, (26) 



d=0 



where £ij(t) and Qj(t) represent the gain effect and the refractory effect, 
respectively. Then, an output of the ijth neuron is defined as follows: 

Xij (t + 1) = /{&(* + 1) + Cij(t + 1)}, (27) 

where f(y) = 1/(1 + e~ v > e ). If Xijit) > 1/2, the ijth neuron fires at time t, 
and the local search to which the neuron corresponds is performed. 

In (|2"5)l . (3 is the positive scaling parameter of the gain effect, and Aij s is the 
gain value of the objective function (|2"4"|) if the local searches are performed. 
Aij S — 9(Sb) ~ 9(Sa), where Sb and Sa are solutions before and after the 
local searches are performed respectively. Here, s indicates a customer to be 
exchanged (Fig. Ufa)) or relocated to their next order (Fig.[D(b)), and s is so 
selected that Ajj S takes the maximum gain. By the gain effect, the neurons 
corresponding to good operations become easy to fire. 

In (|2"B"|) . a is the positive scaling parameter, k r (0 < k r < 1) is the de- 
cay factor, and 9 is the threshold value. Then, the refractory effect inhibits 
firing of the neuron which has just fired; this realizes a memory effect with 
an exponential decay. The strength of the refractory effect gradually decays 
depending on the value of k r . 

In the method of [33H331I3S], the neurons are asynchronously updated. 
When all the neurons are updated, a single iteration is finished. 
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Fig. 4 Construction of chaotic neurons and two local searches used in the chaotic 
search method: (a) Exchange of one customer for another customer and (b) Re- 
location of one customer for another place. If the r/'th neuron fires, (a) the ith 
customer is exchanged (j = 1) for another one customer s or (b) the ith customer 
is relocated ( j — 2) to a next order of customer s 

To evaluate the performance of the proposed method, we solved the 
Solomon's benchmark problems 165] and the Gehring and Homberger bench- 
mark problems [66]. We produced an initial solution using the Braysy con- 
struction heuristic method |67j . Then, we treated time windows as hard 
constraints. The parameters in (|2"5|) - (j2"T)) are set as follows: j3 = 0.04, a — 
0.5, k r = 0.9, 9 = 0.9, and e = 0.02. 

Results are shown in Table |U The simulations are conducted by an Intel 
Core 2 Duo 2GHz computer for 1,000 iterations. Table [2] shows the average 
numbers of vehicles and the average total travel distances (italic) for each 
problem type. These results show that the proposed method provides good 
results. In addition, we compared the proposed method with the other con- 
ventional methods (see [64,43,44,45 for details). The results indicated that 
the proposed method has higher performance than the other conventional 
methods by changing the values of parameter effectively depending on the 
constraints of each problem. 



3.2 Motif Extraction Problems 

To identify an important region in biological sequences, the motif extraction 
problem (MEP) is solved in bioinformatics. In this section, we explain a 
Chaotic Motif Sampler (CMS), that employs chaotic dynamics to solve the 

meps ^msmmmmmhr- 

The definition of the MEP can be mathematically described as follows [7T] : 
we have a biological data set S — {s±, S2, ■•-, %}, where Si is the ith sequence 
(Fig. [5]). Each sequence consists of m, (i = 1,2,..., N) elements. In the case 
of DNA or RNA, the elements correspond to four bases; while in the case of 
protein sequences, they correspond to 20 amino acids. V — {v%,V2, ■ ■ ■ ,vn} 
is a set of motifs, where the length of the motif is L (Fig. [5j). Of course, the 
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Table 2 Results for the 100-customer instances from the Solomon's benchmark 
problems 65 and the other customer instances from the Gehring and Homberger 
benchmark problems [66 . In this table, R means a random allocation, C a clustered 
allocation, and RC their mixture 



type 100 200 400 
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800 
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Rl 12.6 18.4 36.7 
1230.53865.4 9166.1 



55.4 73.0 
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R2 3.1 4.1 8.1 
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11.1 15.1 
13978.5 22639. 



19.2 
3 33250.0 
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840.2 2755.9 7337.2 



58.1 76.7 
14278.7 25310. 



96.1 

8 42376.3 



C2 3.0 6.0 12.1 
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18.2 24.2 
7887.0 12196. 



30.3 
9 17563.1 



RC1 12.0 18.4 
1388.1 3446.5 



36.8 



55.4 73.3 
17453.5 30521. 



90.8 
47038.2 



RC2 3.4 4.7 9.4 
1195.2 2719.4 5654.3 



12.6 16.8 
11519.1 18114. 



19.4 
2 27952.0 



alignment of the motifs and the length of the motifs are unknown. Then, 
the aim of the MEP is to find a set of motifs that maximizes the following 
objective function: 



1 L 



fciuefi 



2 pM ' 



(28) 



where /fe(w) is the appearance frequency of an element (a base in the case of 
DNA or RNA sequences and an amino acid in the case of a protein sequences) 
uj e J? at the fcth position of motif candidates; fl is a set of bases or a set 
of amino acids; and p(u>) is the background probability of appearance of the 
element w (Fig.EJ. In p8|). 01og 2 is denned to be 0. 
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Fig. 5 Definition of the motif extraction problem 



In the CMS, to extract motifs v\, V2, ...,vn, a chaotic neuron is assigned 
to the head positions of all motif candidates (Fig. [6]). The firing of the ijth 
chaotic neuron is then defined as follows: 
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*«(*) = fiVijit)) > ~> (29) 

fW = T^~l M' (30) 

l + exp(-j//ej 

yij(t) = ^(t) + Cij(t), (31) 

where Xij(t) is an output state of the ijth chaotic neuron at time t. If the 
chaotic neuron fires or Xij (t + 1) > h> the zj th motif position becomes a motif 
candidate. On the other hand, if the ijth neuron is resting or Xij(t + 1) < h, 
the ij th motif candidate is not selected, yij (t) is an internal state of the ijth 
chaotic neuron at time t. The internal state of the chaotic neuron |23| is 
decomposed into two parts. The first part fy (t) represents the gain effect of 
the ijth neuron at time t and the second part Qj (t) represents the refractory 
effect of the ijth neuron at time t. They have different effects to determine 
firing of a chaotic neuron in the algorithm. 
The first part &_,■(£) is defined as follows: 

^ ij (t+l)=0(E ij (t)-E), (32) 

^•(') = z££/*m 1o s 2 #t> ( 33 ) 

where j3 (> 0) is the scaling parameter of the gain effect; Eij(i) is the objective 
function when a motif candidate position is moved to the jth position in the 
sequence Sf, E, the entropy score of the current state. If a new motif candidate 
position is better than the current position, the quantity on the right-hand 
side of (f3"2"|) becomes positive; a positive value leads to firing of the neuron. 
Thus, the gain effect encourages firing of the neuron, and such behavior is 
characteristic of a greedy algorithm. 

However, the greedy algorithm gets stuck at a local minimum. To escape 
from the local minima, the refractory effect is assigned to each chaotic neuron. 
The second part Qj(t) realizes the refractory effect. The refractory effect is 
one of the important properties of real biological neurons: once a neuron fires, 
a certain period of time must pass before the neuron can fire again. In the 
model of a chaotic neuron, the second part is expressed as follows: 

t 
Cij (t+l) = -a^ k d r x l3 (t-d) + 9, (34) 

d=0 

= -axij(t) + fcrCtfCO + 0(1 - M, (35) 

where a is the positive parameter; k r , the decay parameter that takes values 
between and 1; and 9, the threshold value. Thus, in (|M|) . Cy(i+1) expresses 
the refractory effect with the factor k r because the more the neuron has fired 
in the past, the more negative is the first term on the right-hand side of (|34[) . 
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Fig. 6 Assignment of the chaotic neurons 

This, in turn, reduces the value of Cij(^ + 1) an d causes the neuron to become 
a relatively resting state. 

Figure [7| shows the results for real biological data set [T!lMHTIl4ini4TJ] . In 
Fig. [7J we show the average probabilities (%) of finding motifs in 50 trials. 
If the motifs are correctly found in 40 trials, the probability is 40/50 = 80%. 
In one trial, we change each motif candidate 500 times. 

From Fig. |7Ja), if j3 takes a small value (J3 = 20), the CMS shows low 
performance. Then, as the value of j3 increase (/3 = 40), the performances of 
the CMS becomes better (Fig.[7Jb)). However, we cannot find motifs for too 
large values of /3 (Fig. |7Jc)). The reason is that the CMS cannot escape from 
local minima if the strength of the greedy effect is stronger than that of the 
refractory effect. In other word, the searching approach is similar to steepest 
descent method. 
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Fig. 7 Results of the chaotic motif sampler 
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3.3 Packet Routing Problems 

A packet routing problem is one of the dynamical combinatorial optimization 
problems because the searching space dynamically changes depending on the 
state of the computer networks. In this section, as one of the applications 
of the chaotic neurodynamics to the dynamical combinatorial optimization 
problems, we explain a packet routing algorithm with chaotic neurodynamics 
for solving the packet routing problems 50, 51, 52, 53, 54, 55, 56, 57, 58 . 



150 



T. Ikeguchi et al. 



Source 



Destination 



Source 




O 



Destination 



Fig. 8 An example of ideal computer networks. In this example, a gray packet 
is transmitted from a source to a destination. An arrow from the source to the 
destination expresses the shortest path of the gray packet. Because there are no 
packet congestion in the ideal computer network, we can easily find the shortest 
path of the gray packet. However, if we simply apply the basic strategy, such as 
the Dijkstra algorithm, to the real computer network, the packet congestion easily 



The packet routing problem is how to transmit the packets to their desti- 
nations as quickly and safely as possible depending on states of the computer 
networks. If a computer network is ideal, the buffer size of each node is infi- 
nite and throughputs of the nodes do not change. In such an ideal case, basic 
algorithms for finding the shortest path length, for example, the Bellman- 
Ford [73], the Dreyfus [74], and the Dijkstra algorithms [75], can find an 
optimal solution of the packet routing problem or the shortest paths for the 
packets. 

However, in the real computer networks, the buffer size is finite and the 
shortest path between any two nodes changes depending on the amount of 
flowing packets in the computer networks or packet congestion. In other 
words, the computer network is one of the dynamic and stochastic net- 
works [7oT7T] . Because the shortest paths between nodes in the dynamic and 
stochastic networks are always changing depending on the state of the net- 
work, we have to consider how to avoid such congestion and how to transmit 
the packets securely and effectively by more sophisticated strategies. 

Now, we define an objective function of the packet routing problem as 
follows; 

r* = min r^ (i = 1, . . . , N g ), (36) 

where N g denotes the number of existing packets; r^ denotes the jth path 
from a source to a destination of the «th packet and it depends on a network 
state; Ri is the set of all possible paths rij. Equation (|3"6"|) means to find r* 
which is the shortest r^ depending on the state of the computer networks. 
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The computer network model has N nodes, and the ith node has Ni adja- 
cent nodes (i = 1, . . . , N). Although there are several ways of how to assign 
a neural network to each node, we take the same way as [T8. 79, 80. 50U55J . 

To realize the packet routing algorithms with chaotic neurodynamics, first, 
we construct a basic neural network which functions to minimize a distance 
of a transmitting packet from the ith node to its destination. To realize this 
routing strategy, an internal state of the ijth neuron, £y , in the basic neural 
network is defined as 

Si i (t + l)=p(l- (kj+d £ {PiW) ), (37) 

where dij is the static distance of a path from the ith node to the jth adjacent 
node; Pi(t) is a transmitted packet from the ith node at the ith iteration; 
g(pi(t)) is a destination of Pi(t); 4?'s(pi(*)) * s ^ ne dynamic distance from the jth 
adjacent node to the destination of pi(t), that is, dj^OpJt)) changes depending 
on Pi(t); d c is the diameter of the computer network; (3 is the normalization 
parameter which takes a positive value. 

If the packets are transmitted to the destinations along only with the 
shortest paths, almost all the packets might be transmitted to the nodes 
through which many shortest paths pass. This behavior might lead to delay 
or lost packets. To avoid such an undesirable situation, one of the possible 
strategies is to memorize a node to which packets have just been transmitted 
for a while, and not to transmit the packets to the node. Then, we use a 
refractory effect peculiar to a chaotic neuron model [H] . The refractory effect 

is defined by 

t 

Cy (t + 1) = -aJ2 k r^ (t-d)+ 9, (38) 

d=0 

where a is the positive control parameter of the refractoriness; k r is the decay 
parameter of the refractoriness and takes between and 1; 9 is the threshold; 
Xij (t) is the output of the ijth neuron at time t which will be defined in (1401) . 
Although the basic mechanism for the memory effect is realized by (|37|) and 
(|38p , mutual connections among neurons are also introduced to control firing 
rates of neurons, because too frequent firing often leads to a fatal situation 
of the packet routing. The internal state of the mutual connection effect is 
described as follows: 

rnj(t + l) = W-Wj2^j(t)^ (39) 

where W is a positive parameter and Nt is the number of adjacent nodes at 
the ith node. Because W > 0, if the number of firing neurons increases, then 
the second term of the right hand side becomes large, which again depresses 
the firing of the neuron at time t + 1 and makes rjij(t + 1) small. 
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Then, the output of the ijth. neuron is finally defined by the sum of the 
above-introduced three internal states, £ij(t + 1), dj(t + 1), and r\ijit + 1) as 
follows: 

Xij (t + 1) = /{& (t + 1) + Cy (t + 1) + %■ (t + 1)}, (40) 

where f(y) = 1/(1 + e^ y ^), and e is a positive but small parameter. In (|i0l . 
if Xy(i + 1) takes the value larger than 1/2, the ijth neuron fires. 

We compared the proposed chaotic routing strategy with a packet routing 
strategy using a neural network which has only the descent downhill dynamics 
of (|37| (the descent downhill routing strategy) and a packet routing strategy 
using a tabu search (the tabu search routing strategy [2"5"ll2"rj| ) . 

We conducted computer simulations of the packet routing for the scale-free 
networks [5Tj. Because real communication networks are scale- free [HI], we 
adopted the scale-free topology as the network topology. 

To evaluate performance of the three routing strategies, we used an arrival 
rate of the packets, and the number of packets arriving at their destinations. 
In this simulation, 20 scale-free networks of 100 nodes are prepared, and 
the quantitative measures, an arrival rate of the packets and the number of 
packets arriving at their destinations, are averaged over these 20 scale-free 
networks. Although we show only the results of the 100-node networks, we 
obtained the similar tendency for other 50-node and 200-node networks. 

Results for the scale-free networks are shown in Fig. [S] In Fig. |H1 the 
proposed chaotic routing strategy keeps the higher arrival rate of the pack- 
ets than those of the descent downhill and the tabu search strategies for 
every packet generating probability (Fig.[9fa)). In addition, the chaotic rout- 
ing strategy transmits more packets to their destinations than the descent 
downhill and the tabu search routing strategies for every packet generating 
probability (Fig.^b)). These results indicate that the chaotic neurodynam- 
ics is effective for avoiding the packet congestion by using the past routing 
history, which is realized by the refractory effect ([58]) . Then, the chaotic rout- 
ing strategy effectively routes the packets to their destinations without loss 
of the packets. 

To reveal the effectiveness of the chaotic neurodynamics for the packet 
routing problems, we analyzed the network behavior of the chaotic routing 
strategy and the descent downhill routing strategy by spatial firing rates of 
neurons. The spatial firing rates by the routing strategies are shown in Fig. 
1 101 which demonstrates that neurons in the chaotic routing strategy (Fig. 
UOH a)) are firing more uniformly than in the descent downhill routing strategy 
(Fig. rrU^b)). Figure ITD1 shows that many paths for the packets are selected 
by the chaotic neurodynamics because the neurons in the chaotic routing 
strategy are uniformly firing during the simulations. As a result, the chaotic 
routing strategy transmits many packets to their destinations by selecting 
the transmission paths for the packets effectively. 

The effectiveness of the chaotic neurodynamics for avoiding congestion of 
packets is analyzed using the method of surrogate data [54.57 . In addition, 
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Fig. 9 Relationship between the packet generating probability and (a) the arrival 
rate of the packets, and (b) the number of packets arriving at their destinations 
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Fig. 10 The spatial firing rates of neurons by (a) the chaotic routing strategy and 
(b) the descent downhill routing strategy. The packet generating probabilities in 
these figures are set to 0.4 

the above-mentioned routing strategy with chaotic neurodynamics is modi- 
fied by adding the waiting times until which the packets are transmitted at 
the adjacent nodes [561158] , The results show that the modified chaotic routing 
strategy exhibits higher performance as compared to the conventional rout- 
ing strategy with chaotic neurodynamics [5U[521[53, 54, 57 and Echenique's 
algorithm [551154] . 



4 Conclusions 



In this chapter, we have reviewed three methods for solving combinatorial 
optimization problems by using chaotic dynamics. 

The first algorithm uses chaotic time series as additive dynamical noise 
that is injected to descent downhill dynamics of the recurrent neural network, 
or the Hopfield-Tank neural network. In this case, chaotic sequences are used 
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to shake internal states of the Hopfield-Tank neural network in order to avoid 
undesirable traps into local minima. 

The second algorithm is also based on the Hopfield-Tank neural networks, 
but the chaotic neuron proposed by Aihara et al. is used as a basic element. 
In this method, the refractoriness produced by the chaotic neuron model is 
effectively used to avoid undesirable local minima. 

The third method used chaotic dynamics to control executions of local 
search algorithm, such as the 2-opt algorithm, the 3-opt algorithm, the Or- 
opt algorithm [19] , the Lin-Kernighan algorithm [20| , and the stem- and- cycle 
ejection chain method p?Tl,n?2] . In this chaotic method [TH[TMi~6U 1 7LfT5] , exe- 
cution of the local search algorithm is encoded by firing of the chaotic neuron. 
Once a chaotic neuron fires, the firing of this chaotic neuron is inhibited for 
a while by the refractoriness, which restricts frequent firing of the chaotic 
neuron, or frequent execution of the same local search algorithm. Thus, the 
chaotic search can escape from local minima efficiently. 

Generally speaking, attractors produced from chaotic dynamical systems 
have fractal structure in the state space, which has the zero Lebesgue mea- 
sure. Thus, effective search using chaotic dynamics can be realized on such 
fractal attractors, which leads to higher performance than those using ran- 
dom dynamics, because the searching space of such fractal attractors are 
much smaller than that of stochastic search [5]. In addition, the algorithms 
using chaotic dynamics can be easily controlled due to its deterministic 
property. 

As we have already shown, the third algorithm in which chaotic dynam- 
ics controls execution of local searches exhibits the best solving performance 
among the three methods. One of the key factors so that the third algo- 
rithm shows the highest performance is that chaotic search is realized with 
the refractory effect, or an exponential decay of the tabu effect. Moreover, 
the algorithm with chaotic dynamics can be easily implemented by analog 
circuits, which can drastically reduce the computational time to obtain good 
solutions. One of the limitations of the chaotic searching methods is that 
we have to tune parameters in the algorithms. However, this drawback can 
be resolved by developing an automatic parameter-tuning method based on 
analyses and controls of chaotic dynamics. 

We have also reviewed applications of the third method to the vehicle 
routing problems [4411431145] , the motif extraction problems [68ll46ll47ll48l l69l 
IHCZD] , and the packet routing problems [5I1IH[51[51[57] . The results of 
computer simulations clearly show that the chaotic dynamics is very effective 
to solve these real- world application problems. 

Although we have only discussed and showed the efficiency of the chaotic 
methods by computer simulations in this chapter, one of the most important 
research directions is to implement these algorithms by analog circuits. By 
the analog-circuit implementation of the chaotic search methods described in 
this chapter, we could develop a novel frontier of information processing that 
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is based on a new computation principle by such nonlinear analog circuits 
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